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Ha es 
ABSTRACT 


The approximation of an ideal frequency response 
by a realizable filter has wide applications in Engineer- 
ing. This topic is treated in this study. The general 
approach used is to convert the problem into the time 
domain and to find a filter satisfying the desired impulse 
response such that the error in the frequency domain is 
minimized. The error criteria used are the Least Integral 
Squared Error and the Maximum Deviation, in the frequency 
domain. The Integral Squared Error (ISE) in frequency 


domain is related to that in time domain by the relation 
oo +00 Oe 
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where F(jw) and H(jw) are the desired and approximate 
frequency responses and f(t) and h(t) are the correspond- 
ing impulse responses respectively. The minimization of 
ISE alone need not limit the maximum deviation in fre- 
quency domain. The deviation in frequency domain is 
related to the average error in time domain by the rela- 


tion 
[F(jw) - H(jw)| < J7?[£(t) - p(t) lat 
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This relation is used to control the frequency domain 
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deviation. It has been shown that the minimization of 
deviation in the frequency domain can be achieved while 
keeping the ISE viehen allowable limits. 

The approximation is carried out using a set of 
orthonormal functions of exponentials. The numerical 
evaluation of the time domain representations of these 
orthonormal functions are carried out by a novel method. 
This method makes the computer evaluation of the transi- 
ent response of any rational function of complex frequency 
simpler and more efficient. 

The theory developed in this dissertation is 
applied to the specific example of approximating the ideal 


low pass filter. 
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CHAPTER 1 


INTRODUCTION 


In recent years there has been a renewed inter- 
est in the approach of approximating a desired frequency 
domain characteristic by converting the problem into one of 
time domain approximation (1-4). This has led to the possi- 
bility of realizing filters using least integral squared 
error criterion in the frequency domain. The family of 
filters having this frequency domain criterion have not 
received wide attention. This is partly due to certain 
problems arising in the solution of the optimum filter 
parameters. This research is intended to fill this gap. 
The general approach employed in this work is to convert 
the problem into a time domain approximation. The least 
Square criterion does not take into account the maximum 
deviation of the achieved response from the desired fre- 
quency response. This work develops a new method of de- 
Sign which minimizes this deviation in the frequency do- 
main while keeping the integral square error within allow- 
able limits. Thus this new method makes use of a combina- 
tion of least square and Chebyshev criterions in the fre- 
quency domain. One problem encountered in this approach 
is the need to compute on digital computer the transient 
response of a Laplace transform expressed as rational func- 


tion. This has been solved by a novel method (5). 
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Having broadly outlined the scope of this work 
in this chapter, the basic problem of least square design 
in the frequency domain is presented in Chapter 2. The 
conversion from frequency domain to time domain and the 
equivalent time domain approximation by using a set of 
orthonormal functions of exponentials are discussed in this 
chapter. The concept of complimentary filter and the 
evaluation of the integral squared error using this fil- 
ter are also reviewed in this chapter. Chapter 3 gives a 
summary of recent works concerned with this problem. 
Chapter 4 gives the mathematical basis of the new method 
of frequency domain design based on minimizing the devia- 
tion in the frequency domain while keeping the integral 
squared error within allowable limits. The basic theory 
of the new method of numerical evaluation of the transi- 
ent response is given in Chapter 5. The convergence pro- 
perties and comparative merits and demerits of this new 
method are also discussed in this chapter. The implemen- 
tation of the new method of design on digital computer 
is done in Chapter 6. The computer algorithm that is used 
to determine the optimum poles of the filter by minimizing 
the integral squared error is reviewed in this chapter. 
The results obtained by applying this new technique to a 
specific example are also discussed here. Chapter 7 gives 


a summary of the entire dissertation and discusses the 
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conclusions of this research work. The areas of further 


work are also mentioned in this chapter. 


CHAPTER 2 


LEAST INTEGRAL SQUARED ERROR FILTERS 
AND 


EVALUATION OF INTEGRAL SQUARED ERROR 


INTRODUCTION 

In this chapter we define "Least Integral Squared 
Error Filters" and outline the technique of realizing such 
filters. The theory of evaluating the "Integral Squared 
Error" (ISE) by a simple filtering operation is also dis- 
cussed, 

Let us consider the problem of approximating an 
ideal frequency response F(jw) by a physically realizable 
filter. The impulse response of any physically realizable 
filter can exist only for positive values of time (6). Let 
H(jw) be the frequency response of such a filter. If f(t) 
and h(t) are the inverse Fourier Transforms of F(jw) and 


H(jw) respectively, we have 
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F(jw) = Sp £(t)e 2? Fat 
+00 4 
and H(ju) = J h(t)e I°*at 
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Let E(jw) be the error in the frequency domain and e(t) 
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be the error in time domain. 


Then 


Hi) eB Cte 00) 


and e (t) f(t) - h(t) (222) 


There exists a vast body of literature on approximating 
specific examples of F(jw) by making use of purely fre- 
quency domain techniques (7-10). This research develops 
a general technique of design of filters that can be 

adopted to any specific example without difficulty. It 
May also find application in the design of optimum fil- 


ters in signal theory (11). 


LEAST INTEGRAL SQUARED ERROR FILTERS 

All methods of approximating F(jw) by H(jw) try 
to minimize the error E(jw) in equation (2.2) in some 
sense, The family of filters having minimum integral 
Squared error in frequency domain has not received atten- 
tion until recently (1-2). This work develops a method of 
design based on minimizing the integral squared error 
f_ |B (ju) |2au. The filters satisfying this criterion are 
defined as Least Integral Squared Error Filters. 

The integral squared error (ISE) in frequency 


domain and time domain are related by Parsevel's theorem 
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(12, 16). It can be easily shown that, if e(t) belongs 


tO, E(-o, 0), 


+a eo 
Sof |E(jw)|7dw = f*ele(t) [Pat OR) 
where e(t) and E(jw) are defined in equation (2-2). The 


function be) ein equation (2,2) exists only for t =. 0. 
Hence it cannot be used to approximate f(t) for t < 0. 
Therefore it is assumed that f(t) vanishes for negative 


values of t. Equation (2.3) can now be written as 
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acl |E(ju) | 2aw = for je(t) |Zat (a4) 


Equation (2.4) converts the problem of least square design 
in the frequency domain to a problem of least square de- 
Sign in time domain. This makes it possible to use time 
domain approximation theory for frequency domain design. 
If £(t), the inverse of F(jw), is known the problem of 
design in the least ISE sense reduces to finding h(t) 
which minimizes the expression fo’ je(t) |Zat. Mies cRroOLce 
of h(t) is governed by the factor that the filter is to 
be realized from passive, lumped linear circuit elements. 
This requires that H(s), the Laplace transform of h(t), 
be a rational function of s with all the poles lying in 


the left half s plane. The impulse response of such cir- 
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cuits can only be a finite sum of damped exponentials. 
(The case of multiple poles is exceptional and is not con- 
sidered.) This suggests that the best form of approxima- 


tion for f(t) is by a sum of damped exponentials. 


EXPONENTIAL APPROXIMATION 
Let £(t) be approximated by h(t) which is a lin- 
ear combination of N damped exponentials. 
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h(t) = fee a= a8 C29 


Ss; is, in general, a complex number with Re (s,)<0. Let I 
be the time domain integral squared error. 
Then 

t= tp (etc 2 ayei*y7ae (2.6) 


i= 


The parameters (ee and (sens have to be determined 


=1 
such that I is a minimum. This problem was first dis- 
cussed by Aigrain and Williams and they formulated a set 
of equations known as Aigrain and Williams equations (13). 
These equations are of classical interest and are neces- 
sary to understand the complexity involved in exponential 


approximation. 


Differentiating the equation (2.6) with respect 
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For I to be minimum, 
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It is possible to simplify the equations (2.7-2.9) as 
follows. 


By definition, 
ee (tye *k*at = F(-s,) (2.10) 


where F(s) is the Laplace transform of f(t). 
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and 


- i 
dt = + 5," ao 


Substituting equations (2.10) and (2.11) in equations 


(2.7) and (2.9) we get, 


N ee 
Ea — — _— 
ret (enriesey = Bi Sy) k = Lie prarece (214) 


Similarly substituting equations (2.12) and (2.13) in 


equations ((2. 8) “and (2.9), 


= 8 
a poe = -F'(-s,), k = 1,2,...N Carl) 
i=l \°i k 


Equations (2.14) and (2.15) are known as Aigrain and 
Williams equations (13). There are 2N unknown complex 
parameters and 2N equations. The analytical solution of 
these 2N nonlinear simultaneous equations is extremely 
difficult; but these equations suggest the probable exis- 
tence of a set of parameters co ae and (ena which 
will make the ISE, I, a minimum. Moreover any direct 
solution of these equations requires the availability of 


F(s) and F'(s) at any point. Because of these factors no 
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attempt is made to solve these equations. Young and 
Huggins have shown that the complexity involved in tack- 
ling these equations could be greatly reduced by making 


use of orthonormal functions of exponentials (18-21). 


ORTHONORMAL FUNCTIONS OF EXPONENTIALS 


Let to, (t) } be a set of orthonormal functions 


i=l 


defined over (0, ~). Then 


Foo, (t) >, (t)dt 1,.if i= j 
(2.16) 


0, if i # 3 


Any veal function £(t), te(0,~), for which | £(t) |Zat < © 
can be expanded in terms of these orthonormal functions. 
If the set of orthonormal functions is complete and if 

the series E0405 (8) converges uniformly, we have (12, 


i= 
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where 
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h£(t), (t)dt (2218) 


If f(t) is approximated by N terms, the approximation h(t) 


is given by 


h(t) = q COC tc) (2.19) 
i=l 
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The error of approximation is 


e(t) F(t)erehte) 


N 
E(t) - 25C44; (t) (2.20) 


It can be shown that ISE is given by, (14), 
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It is possible to construct a set of orthonormal functions 
of exponentials by Gram-Schmidt process (12, 14). Another 
method of greater practical value was developed by Kautz 
(17) which was later improved by Young and Huggins (18- 
20) and Ross (21). This has been widely used in the ana- 
lysis of speech signals (22-23) but has not been adapted 
and made use of in filter design. This dissertation makes 
use of these orthonormal functions. 

Let it be required to construct a set of N ortho- 
normal functions from N exponentials of which there are 
r real exponentials and c pairs of complex conjugate ex- 
ponentials. It is assumed that all the exponentials have 
negative real parts. 

Consider the frequency domain representations 


of N functions as given below. 
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where 


fa 


q.= |s,_,|7 = |s 
Cc N-l N 


Sy-l and Sy form the last pair of complex conjugate expo- 


nentials. 


Functions oF to - correspond to the set of real exponen- 


tials, eerie and 9-42 to the first pair of complex exponen- 

tials and PHL and o to the last pair of complex exponen- 

tials. It will now be proved that the corresponding time 
N 


domain functions to, (t)} form an orthonormal set. 


i=l 


This is done by using Parsevel's theorem. By this theo- 


rem we have, 
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Consider two functions #4 
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corresponding to a pair of complex conjugate exponentials 


ne eel and Sr42i° Using equation (2.22) we get, 
2p,vq, s 
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The evaluation of yb 2Se (s)¢@ 


yore oieT (-s)ds may be done 


rel 
by taking the closed contour consisting of the imaginary 
axis and left side semicircle at infinity. Thus by 


Cauchy residue theorem (25) we have 


ie Dant (s)@ 


Sse i eee ee ee em Oar 26) 
where R, and R, are the residues of {%,,,_)(S)9,5, (-s)} 
she ee ee and S42 respectively. Evaluating the resi- 


dues it can be easily shown that (Ry + R,) is zero. This 


proves that functions $9 (Cojleand @ (t) are orthog- 
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onal to each other. Similarly 
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Integrating these functions along the closed contour enclo- 
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Thus the orthonormality of the functions d4o4-1] (t) and 
> 424 (t) is proved, 

If any two functions bats) and o,(s), which do 
not belong to the special case discussed above, are con- 
sidered it can be shown by pole-zero cancellation, that 
all the poles of 9, (8) ee) lie only on one half (left 
or right) of s-plane and the function is hence analytic 
in the opposite half of s-plane. Therefore it can be 
proved that 


paleo = me 
ie een on s)ds = 0 (230) 


Considering any function o.(s), which is formed from real 


poles only, we have 
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®, (Ss) ®, (-s) = Daa} (2.381) 


Evaluating the countour integral by Cauchy Residue theorem 


we Can prove 
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Thus the orthonormality of functions ee) eee is proved. 
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COMPLIMENTARY FILTER AND ISE 
When ISE is minimum with respect to the para- 


meters {8., a, 3 of “equation (2.5) 
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and lee 0 
ds. 
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From equations (2.7) and (2.8) we get 
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Equation (2.33) gives the condition for I to be a minimum 
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In other words, the condition for I to be minimum is that 


e(t) should be orthogonal to tone nm and Gee wore 


These properties of e(t) are made use of by Young and 


Huggins (18-20) in developing a simple method of evaluat- 


ing ISE. This is now discussed. 
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In the formation of the orthonormal functions 
o(s), it is observed that each %(s) has an all pass 
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tor successively increases, If N exponentials are used 
for the approximation and it is required to add one real 
pole or a pair of complex poles, the all pass filter 
appearing in the corresponding orthonormal function or 
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called the complimentary filter of Nth degree approxima- 
tion (18-20). This filter has the interesting property 
that if the function f(t) is time reversed and filtered 
using this complimentary filter, say G(s), the integrated 
Squared error can be directly evaluated from the output 


of G(s) (18-20). This property can be proved as follows. 
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Since h(t) can also be expressed as a sum of exponentials 


as in equation (2'.5)), 
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h(-t) = 2 B.e (2.38) 


Comparing equations (2.37) and (2.38) it is seen that the 
approximation of v(t) can be done solely on negative time 
using growing exponentials as against the decaying expon- 
entials used for approximation in positive time. In the 
frequency domain the corresponding functions will be 

{, (-s)}4_) instead of fas) i, 2s Two-sided Laplace 
transform is used to get the frequency domain functions 
of negative time functions (24). 
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two terms A'(s) and A''(s) such that 
A(s) = A'(s) + A'"(s) (2.50) 


where A'(s) = G(s)€(s) has all poles in the right half s- 
plane and A''(s) = G(s) 25,9; (s) has all poles in the 
left half s-plane. The output a(t), of the complimentary 
filter, with reversed f(t) as input, extends from -~ to 
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a'(t) and a''(t) correspond to A'(s) and A''(s) respec- 
tively. It is easily seen that a'(t) is the output of 
the complimentary filter during negative time and a'' (t) 
during positive time. 
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where a(t) is the output of the complimentary filter with 


reversed f(t) as input. 


This is schematically represented in Fig. (2.1). Methods 
of implementing this on digital computer are discussed in 


Chapter “6r, 
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Figure 2.1. Schematic diagram for evaluating ISE 
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N parameters, instead of the original 2N parameters of 


equations (2.14) and (2.15). 
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CHAPTER 3 
A SURVEY OF PREVIOUS RESEARCH 


In this chapter we discuss some of the impor- 
tant methods that have been used to determine the para- 
meters of H(s) in the least square sense. Some of these 
methods tackle the problem purely as a time domain approx- 
imation, while others take the frequency domain behaviour 


into consideration. 


PRONY'S METHOD 
The first attempt to approximate a time function 
by the sum of weighted exponentials was due to Prony (26). 


Prony chose the function 
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The basic philosophy of Prony's approach is made 
clear by equations (3.9) and (3.5). Equation (3.9) defines 
a difference equation and (3.5) corresponds to the charac- 
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the coefficients of an equivalent difference or differen- 
tial equation is common to most of the methods reviewed 

in this chapter. Prony's method has also been extended to 
the case where the available number of points of f(t) is 


more than 2N, the number of unknown parameters. 


KAUTZ'S METHOD (17) 

tf £(t), the function to be approximated is the 
sum of N weighted exponentials, then f(t) is the solution 
of an Nth order, constant coefficient linear differential 


equation. This means 
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Here e(t).<is not the exact error of approximation, but is 
a measure of how closely f(t) can be expressed as a sum 


of N exponentials. The eye 
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these points of discontinuity are undefined and hence the 
method has to be modified. Another point worth noting is 
that, instead of a = 1 in equation (3.12) any of the a. 
can be taken as unity and this will lead to a different 
set of solutions for the exponentials (26). An interest- 
ing review of this and similar methods, which basically 


make use of the approach due to Prony, can be found in 


McDonough (26). 


YENGST'S METHOD (27) 

Yengst's method resembles the approach of Prony 
in that he assumes a difference equation, the solution of 
which gives h(t), which approximates the desired function 
f(t). The Nth order difference equation used by Yengst 
is 
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exponentials. If the matrix is non singular, equation 


(3.19) can be solved for Ga 


Gel The exponentials are 


then obtained from the solution of the polynomial equa- 


tion 
12 Se eG ay_12 + a = 0 (rey 
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Yengst (27) has shown that (e 


are the solutions of 


Mie] 
the equation (3.20). It may be noted that equation (3.20) 
corresponds to the characteristic equation of the differ- 
ence equation (3.15). 

In Yengst's method, the zeros of H(s) are deter- 
mined by the successive application of the initial value 
theorem to obtain expressions of initial values of h(t) 
and its (N-l1) derivatives at origin. These are then 
equated to the corresponding values of f(t). The solution 
of the N simultaneous linear equations, thus formed, gives 
the zeros of H(s). This method is easy to apply if f(t) 
is known analytically. If £(t) is specified as samples, 
Yengst suggests a polynomial interpolation of f(t) at 


t = 0. This method of determination of zeros is not very 


convenient when f(t) is not known analytically. 


METHOD OF PERDIKARIS AND LAGO 
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time domain synthesis, is very closely related to the 
method of Yengst. This paper makes use of Z-domain 
approach instead of the s-domain approach of Yengst. An 


Nth order difference equation is first assumed. 
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Equation (3.21) can be expressed as 
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By making use of the properties of Z-transforms (29) we 
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where z(y.) = Y(Z). 
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If we consider H(Z) as the equivalent Z-transform of H(s) 


we have 
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The error criterion used by Perdikaris and Lago for deter- 
mining Ten is exactly the same as used by Yengst. 
N 
Hence the error defined as (yoo leary. \e leads to 
i=N 1 kel & =k 
the same set of equations as those arrived at by Yengst in 
equations (3019). “Once (apinet are known the poles in the 


Z-domain are determined and the corresponding poles in s- 


domain are found from the relation 


ak 
Se = Plog Zs (3.26) 


The Z-transform expression of equation (3.25) is made use 


of to find the zeros of H(s). Since we have the relation 
H(Z) > H(s) CSR 29) 


and H(Z) is completely known from equation (3.25), the 

zeros of H(s) can be found by splitting H(Z) into partial 
fractions. The typical factors appearing in the partial 
fraction expansion of H(Z) and their s-domain equivalents 


are as below. 
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By making use of the partial fraction expansion of H(Z) 
and equation (3.28), H(s) can be found. 

An alternate method of determining zeros of 
H(s), suggested by Perdikaris and Lago, is to force h(t) 
to pass through N points as specified by the desired 
function f(t). Since the poles of H(s) are already known, 
h(t) can be expressed as a sum of damped sine, cosine and 
exponential terms. The coefficients of these terms con- 
stitute the N unknowns which determine the zeros of H(s). 


Hence N equations are formed by choosing h(t) such that 
h(t, ) = f(t.), ee Ly 25 oan ra 629) 


The solution of equation (3.29) gives these coefficients 
and the zeros of H(s) are found from them. This method 
has the advantage that h(t) can be made to approximate 


f(t) very closely in a short range but h(t) may deviate 
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METHOD OF VASILU 

A different approach of time domain synthesis is 
found in a paper by Vasilu (30). In this method the im- 
pulse response of a pulse transfer function H(Z) approxi- 
mates the desired function £(t). The poles of the func- 
tion H(Z) are assumed to be uniformly distributed on a 
circle of radius smaller than unity. Thus H(Z) has the 


form 


4(2Z) = (3530) 


where |a)| < l. 


The choice of ap is made arbitrarily. It is proved in 
(30) that PESE (tb) Mis zerosforvte> To ,rtheserrorsin, the 
range t > Ty can be reduced by choosing ap very small. 
But too small a value for a, will result in larger errors 
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the mean square error between the impulse response of 
H(Z) and the desired function f(t) at sampling instants 


is minimized. 
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METHODS OF McDONOUGH AND HUGGINS (26-31) 

All the methods discussed so far do not minimize 
the actual ISE as defined by equation (2.6). We will now 
review some of the recent methods which minimize the ISE. 

McDonough and Huggins make use of the complimen- 
tary filter to find a set of optimum poles. From equation 


(2.52) we have 
te 2 es ai 2 
Deeg (t= h(t) Sates eet). ot (2231) 


where a(t) is the output of the complimentary filter G(s) 


when the reversed signal f(-t) is applied as an input. 
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Figure 3.1 Evaluation of ISE 


Figure 3.1 shows how to evaluate I for a set of 
N exponentials. One method of computing the exponentials 
by minimizing I is to use gradient techniques. But, it 
is found that error I is very insensitive to pole positions 
over a wide range around the optimum poles (26). The 
values of the gradients at these points are nearly zeros. 
Hence gradient techniques are not suited to minimize I. 
This has led to the development of various iteration 


schemes of computing the best pole positions. These are 
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now discussed. 


The condition for I to be minimum with respect 


to (emer and (on are expressed as (equation (2.34)) 
rpe(tjesi*at = 0 (2e2a) 
ree (t)te®itat = 0 (3.33) 

Re ype pec 


McDonough and Huggins (26, 31) make use of equations (3.32) 
and (3.33) to develop an iterative scheme as follows 


Let £(t), defined over (0, ~) be a signal vector 


SaUe yy 
i=1 
subspace Sy Ofec.2) Hquaction. (3.32) simpisessthat Lomsk to 


be minimum the projection of e(t) on Sy should be the null 


in a Hilbert space S. The exponentials (e Span a 


vector. Let us now consider a subspace S of Sesuch that 


2N 
Ss-t,N S-t,N 


S is spanned by (e 1 ee and (te 1 ree Equations 


2N 
(332),-ands(3.33J. ply that, tor. Co .be mMinimum,. e() 
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Hence we have, 


Projection of e(t) on Son - Sy 
= (Projection of f(t) on Son - Sy) = (Projection, of 
levee teye Son - Sy) 


Since h(t) lies entirely on S h(t) has no com- 


N?’ 


ponent on S - § Thus we get a modified condition for 


2N N* 
optimality of the poles. This condition is that, for ISE 


to be minimum, f(t) has to be orthogonal to the subspace 


McDonough and Huggins chose a set of orthonormal 


functions as bases for space Son ~ Sy as follows (31). 
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and ®,(s), i = 1,2,...N are as defined by equations 2.22. 
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Hence o, (t), Kes sl) 2, one NeLormt an ortnonormal 
set and d. (t), k = N+l, ..., 2N are bases for the sub- 


space Son = Sue 


The condition for £(t) to be orthogonal to the subspace 


Son Sy can now be stated as 


Sof (t) > (t)dt = 0, i = 1,2,...N (3.36) 
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where P44 (8) = G(s)?, (s) 


Equation (3.36) can be written as follows 
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The right hand side of equation (3.39) can be described as 


a filtering operation as shown in Figure 3.2 
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Since P44 (8) = G(s), (s), the filtering operation of 


Fig. 3.2 can be modified as shown in Fig. (3.3). 
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v(t)=£ (-t) Say) a Deora cn 
ae re, 
Figure 3.3 Simplified diagram to evaluate ort) $. 2 (trae 
by filtering operation 
The optimality condition is stated as 
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McDonough and Huggins later modified their method 


into an iterative scheme (31). 
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where F(s) is the Laplace transform of f(t). 
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into D(s) and o(s), Ko] li 2yeecyNa  EQuUatLOn (444) scan 
now be solved for Sas 1. = 0,1,7...(N-1)= When equation 
(3.45) is satisfied the iteration is stopped. Otherwise 
these S5 values are used to make a new estimate of (Sryaae 
and a new set of values of g, are found and the iteration 


is continued. 


The evaluation of expressions like 
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can be done by the filtering scheme similar to that in 


Fig. 3.2. Some numerical difficulties are reported in the 
st 
D(s) 


The initial estimate of pole positions has to be near the 


realization of filters and (Ss) in cascade (31). 
optimum to achieve convergence. When f(t) is nearly 


exponential the convergence is fast. 


SEAR'S METHOD 


An iterative scheme, similar to the one dis- 
cussed above, was first suggested by Sears (32). 
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Equation (3.50) may be expressed as 
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Or 
Nat SS) oe ae = 
See at_.b 1¥; (s)F(-s)}L"{G, (s)F(-s)} dt 
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values of ass i = 0,1,...,N-l1 are at the optimum point, 
equation (3.55) is satisfied. If not, the solution of 
this equation gives a new estimate of Coane 
ation of terms like s_ ob 1{P (-s)¥, (s) 2" {F(-s)G, (s) }at 


The evalu- 


is carried out as follows. 


Let 
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(3.57) 
and L(G, (s)F(-s)} = g,¢(t) 


V5 _(t) and Sy. (t) can be generated by the filtering oper- 
ations shown in Fig. 3.4. The evaluation of the expres- 


sion is done by integrating the product of bi ¢ (t) and Sig l(t). 
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Figure 3.4. Schematic diagram to evaluate v; ¢ (t) and 
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the iteration is repeated. Hence the success of the 
method depends on making a 'good' initial estimate of the 
pole positions. 

The methods of McDonough, Huggins and Sears are 
very important. None of the previous techniques reported 
in the literature sought to find an exact solution of the 
exponential representation of signals by minimizing the 
actual ISE. The new methods are computationally more 
tedious than the previous methods, but this is inherent 
in the very nature of the problem as indicated by Aigrain 
and Williams equations. 

In all of these new methods the zeros of H(s) 
are found by first determining the coefficients of the 
orthonormal expansion, as given by equation (2.18), that 
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Ce = fof(t)o, (t)dt 


The evaluation of C. is done by a filtering operation sim- 
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Figure 3.5. Evaluation of orthonormal coefficients 
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Once C. values are found, 


N 
H(s) = 2 C2. o.x(s) (3558) 
Equation (3.58) gives the zeros of H(s). 


MSS METHOD 

An alternate approach to arrive at an exact 
solution was suggested by McBridge, Schaefgen and Steiglitz 
in (34). (For convenience this will be called MSS method). 


In this method, a rational Laplace transform H(s) is first 
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The ISE is given by 
I= sol£(t) - n(t)|7at (3.60) 


where h(t) is the impulse response of H(s). 
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An analytical solution of equation (3.64) is difficult 
because of the nonlinearity. An iterative solution was 
suggested in MSS method by defining a new error e, (t) 

such that e, (t) tends to e(t) as the final solution is 


approached. The new error e, (t) is defined as 
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E.(s) = Bayt 2) - see (3.65) 
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where j is the number of iteration. When final solution 
is reached ES = 2s ak) and E, (s) = Hts). A newrlor 
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In deriving equation (3.67) it is assumed that the para- 
h 


meters of Pha ee are not varied in sthe 5° iteration. 
E, (s) can now be written as, 
eee Sa eae eon pest bse ee 
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and Be foPag (t)p(t)dt Vo ae Ne Vector, (30 13: 
Equating Grad q = 0 in equation (3.72) an iteration equa- 
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Pag (t) is also replaced by eles Corresponding changes 
are made in the matrix P and vector C. 

The second method of MSS cannot be started with 
an initial estimate of zero or unity for H(s). In general 
it is found that method 1 is well suited for initial points 
far from optimum and method 2 converges faster near the 
optimum. It should be pointed out that in MSS method the 
iteration is to be carried out for 2N variables while in 


previous methods only N variables needed to be considered. 


MILLER'S METHOD (35-36) 

Miller modifies the MSS method by introducing 
Aigrain and Williams equations to the error E, (s). Using 
equations (2.7), (2.8) and (2.9) we obtain a different 


form for Aigrain and Williams equations. These are 
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where (saae give the exponentials for minimum ISE. Sub- 
stituting E, (s) and E, ' (s) from equation (3.65) into equa- 
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os (-s,)F(-s, ) = N, (-s,) 
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Baa) used in equation (3.65) does not appear in equa- 
tvon (3.73). Miller has shown (35) that equation (3.7/8) 
can be directly derived from Aigrain and Williams equa- 
cLOnS as Givens by equation, (3.76). ih Lares are con- 
Sidered as an estimate of the pole positions, equations 


(3.78) gives 2N simultaneous linear equations in 2N un- 


knowns (a,, siete The solution of equation (3.78) thus 
gives a new estimate of the pole positions. Iteration is 


continued until there is no further change in pole posi- 
TlLons. 

Miller's approach could be used only when F(s) 
and F'(s) are known. When f(t) is analytically known 
this method is useful. Miller establishes the link between 
MSS method and Aigrain and Williams equations. Equations 
(3.78) gives an iterative method of solving the Aigrain 


and Williams equations. Both MSS and Miller methods con- 
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verge slowly (36). 


APPROXIMATION TO IDEAL LOW PASS FILTER USING TIME DOMAIN 
APPROACH 

Two different methods of approximating the ideal 
low-pass filter, by time domain techniques have been 


reported in the literature (1-3). In (3) Ulstad approxi- 


Since 
15 


of damped sine and cosine functions. The techniques of 


mates the delayed and truncated function by the sum 
nonlinear programming are used to find the sine and co- 
Sine functions and their weights such that ISE is mini- 
mized. These values are then used as the starting point 
for minimizing the Chebyshev error in time domain. 

In a recent paper, Pottle and Wong make use of 
time domain techniques to achieve optimum least-squares 
approximations to ideal low-pass filter (2). The approxi- 
mation is by means of orthonormal functions dy, (t), k=l, 
2,-..N discussed in Chapter 2. The ISE in the frequency 


domain and time domain are given as follows 


| 
| 


1 .+0 : ; 2 
= Suan cd = H (jw) | dw 


pt? | £(e) - w(t) (“dt (3.79) 


F(jw) is the desired frequency response and H(jw) is the 
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approximation. The functions f(t) and h(t) are the in- 


verse of F(jw) and H(jw) respectively. 


and 


The 


The 


(2) 


The 


N 
h(t) = © C.9, (t) 
i=1 
(3°3:80)) 
H(jw) = 2,9, (ju) 
jw) = -®. (Jw 
ethene 
C, = Sof (t)o,; (t)dt 
= 5/7 °F (-jw)o, (jw) (3.81) 
Spree jw); (jw)da ° 
expression for I can be written as 
©0 N 
I= a |F (jw) |7aw - Ge (35:82) 
ip eed Aten tL: 
i=1 
frequency response considered by Pottle and Wong in 
is 
F(jw) =e 2°, fol < 1 
(335) 
= 0, Weigle 23 ok 
inverse £(t) of F(jw) is 
sin (te. =<) 
peal 0 
f(b) = Te (3.84) 


Equations (3.81) and (3.82) can be written as 
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Tot Laie : 
C; = eae 09, (jw) dw (SoD) 
N 
and me ee (3.86) 
T : st 
1=1 


The problem of choosing the parameters of H(s) by minimiz- 
ing I in equation (3.86) was first discussed by Pottle and 
Thorp in (1) and recently by Pottle and Wong in (2). Pottle 
and Thorp make use of the steepest descent technique and 
Newton's method of minimization. The optimization is done 
with respect to the parameters Ps and qs of functions 


o. (jo), i=, Bee We SPhetdélay*t.\ in equation -(se83) 


0 
may be kept constant or may be considered as another para- 
mMeEcver, 

The results obtained by Pottle and Thorp by 
applying this minimization: technique are not Satistactory. 
It is found that in most cases reported in (1) the optimum 
point has multiple poles. This, obviously, is a contradic- 
tion because the orthonormal functions o,(t), ee ree NN. 
Sake 

i 


are only linear combinations of (e and hence , (s) 


=1 
1 =) ee, Peep cannot “have mui biplespoles., 
An improved version of this method appears in 
the recent paper by Pottle and Wong (2). The general 
approach in (2) is the same as in (1). The minimization 


is carried out by using Fletcher-Powell algorithm. The 


optimum points, as reported in (2), do not have multiple 
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poles. 


The following points may be noted in connection 


with the methods in (1) and (2). 


The function f(t) in equation (3.79) exists for -@ < 
t < », The orthonormal functions o;(t), nl eR Pane Cores 
N and h(t) exist only for t 2 0. Hence the condition 
of completeness, required for any orthonormal expan- 


sion, is not satisfied (14). 


The ISE evaluated by equation (3.82) corresponds to 


the total error for -© < ‘t < ©, 


The coefficients ceil ey correspond only to that 


Dale wOfet (talon which to 0. 


Whenever a pair of poles occur near the imaginary 
axis, ®, (jw) will have peaks at the corresponding 
values of w. This leads to serious difficulties in 
the numerical evaluation of (ca and the gradients 


ec) te 


The paper by Pottle and Thorp (1) is the first 


attempt, reported in the literature, to approximate a 


desired frequency response in the least integral squared 
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error sense. This method and its later modification by 
Pottle and Wong (2) have the advantage that when the sig- 
nal is band limited, all integrations over the semi infin- 
ite time axis could be converted to finite integrations 

in the frequency domain. But the numerical difficulties 
reported in both” (il) and= (2)"limit the applicatrons “of 

the method. 

All the methods reviewed in this chapter approx- 
imate a given impulse response by the inverse of a rational 
Laplace transform H(s). In most cases the poles and 
zeros of H(s) are determined in the least square sense. 

In some cases, the poles are located first and the zeros 
are then chosen to satisfy certain conditions in the time 
domain. The determination of optimum poles always involve 
tedious computations which sometimes lead to numerical 
difficulties. The problem of finding H(s) so as to sat- 
isfy specific conditions in the frequency domain, in addi- 
tion to the conditions on ISE, has not yet been reported 


in the literature. 
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CHAPTER 4 


IMPROVING THE FREQUENCY RESPONSE BY MINIMIZING 


AVERAGE ERROR IN TIME DOMAIN 


It has already been shown that minimizing the 
ISE in the time domain is equivalent to minimizing the 
ISE in the frequency domain. From the point of view of 
the frequency domain performance the minimization of ISE 
alone need not produce the best result. For example, at 
some value w, the magnitude response of H(jw) may have a 
sharp deviation from the desired response even if the ISE 
is a minimum (2). This chapter develops a new method of 
Minimizing this deviation, keeping the ISE within allow- 


able limits. 


PROBLEM STATEMENT 
Let F(jw) be the ideal frequency response to be 
approximated and H(jw) be the approximation. The devia- 


tion E(jw) for -™ < Ww < ~ is defined as 
E(jw) = F(jw) - H(jw) (At) 


The problem is to find the parameters of H(jw) such that 
the upper bound of |E(jw)|, we(-~, ~) is minimized sub- 
ject to the condition that the ISE is less than a pre- 


assigned value. The problem as stated above considers 
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both the ISE and the deviations in the frequency response. 
The solution of this problem, thus offers a compromise 
between the least square and Chebyshev criteria in the 


frequency domain. 


AN UPPER BOUND ON DEVIATION 


The error E(jw) is defined as 


E (jw) FG) 0) 


where e (t) 2a Gad oc lel Ged) (4.2) 


E(jw) can be expressed as 


E(jw) = Ste (t)e J°tat 


—=0oO 
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Since e(t) eL*, we can apply Schwarz's inequality giving 
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Thus we have 
JE(jw)| < S72 le(t)|dt (4.3) 


The function h(t) of equation (4.2) is formed from a set 


of N exponentials as discussed in Chapter 2. Hence 


N 
h(t) = 2 Cz, (t) (4.4) 


where >, (t), 1 = 1,2,...,N are the N orthonormal functions. 
Using equation (4.4) equation (4.2) can be expressed as 
N 


Ct) ea Cc) C.o, (t) (4.5) 


i=l 
Therefore 
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Each of the functions o.(t), pele 2 pare See SUMO 
damped exponentials and exists only for t 2 0. From Chap- 


ter 2, we have 
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6, (t) 


| 

I ms 
ve) 
0) 


{Re(s,) + jIm(s.)}t 
R.e = ; ‘ (4.7) 


where R, is a complex constant andl<k<«<«n<«<N. 


ae eae (Re(s2) (oj im(s.) te 
1L= 
n 
Sagi (Ry [yoes> 91) tae 
i=l 
e il 
~ 2, a! =e TEST ae 


By using equation (4.8) we get 
N J 
Z (Ci [fo], (t) [dt < (4.9) 

i=l 


(assuming each C; to be finite) 


Using equations (4.6) and (4.9) we have the result that 

fi lett) [at is finite, when pr fe (t) lat is finite. This 
leads to the conclusion that the deviation E(jw) has a 
finite upper bound when fy cei tat exists. Moreover if 
e(t) is unidirectional fe ete wae is the same as the mag- 


nitude of the deviation at w = 0. Hence, in the most gen- 


eral case sr? le (t) [dt is the least upper bound of |E(jw) 


This expression, fe betty aay is,derined as Che average 
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error of approximation. Thus the minimization of average 
error offers a method of minimizing the deviation in the 
frequency domain. This is now considered in greater 
detail. 

Since >; (t) = 0 onset y= 0 equation: (4. @)can: be 


modified as follows. 


srBle(t) [at = sTBle(t) - 2 c.g, (t) fat 


| end 
ltrs 
= 


y_|£(t) Jat + folf(t) - . C$, (t)|dt (4.10) 
i=l 
Examining equation (4.10) we see that the choice 
Of Cio; (t), i = 1,2,...N has no influence over the term 
f Ve (e) fat. For a Causal system £(t) ==) 0 for Ea 0) and 
this term vanishes. Hence in order to minimize the aver- 
se error we need only consider the expression ARTES: = 


5 C.$, (t) | dt. Let y represent this average error 
i=l ~ 


N 
Yh = A at| HGy oe SOMO ME Fae (4,11) 
nal 


The problem of minimizing the upper-bound of 
deviation in frequency domain subject to a constraint on 
ISE can now be restated as follows. 


Choose the parameters of h(t 
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minimized subject to the condition that the integral squared 
N 

E6rrorel a= fol f(t) = ae Cio, (t)|7at < U, where U is the 
Til 


allowable ISE. 


MINIMIZATION OF AVERAGE ERROR 

The determination of the parameters of H(s) = 
10a 8; (8) is carried out in two steps. The first step 
is to find a set of poles of H(s) which will minimize the 
ISE. These poles are the solution of the Aigrain and 
Williams equations (2.14-2.15). Even though the unique- 
ness of the solutions of equations (2.14-2.15) has not 
yet been proved mathematically, all the researches re- 
viewed in Chapter 3 have reported that in all of the cases 
considered the various iterative schemes have converged 
to the same point. Hence it may be assumed that the solu- 
tion of the Aigrain and Williams equations or equivalent 
equations gives a set of poles yielding the least ISE. 
Once the poles ee LN are known a set of orthonormal 


functions >; (t), 1 = 1,2),..s.N can be formed. The second 


N 


step is the determination of the coefficients (Ci) 56] 


such that y is minimized, subject to the constraint on 
TSE. 
Let Epes be the set of poles giving the least 


ISE and >, (t), i= 1,2,...N be the corresponding ortho- 
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defined by equation (4.22), forms a compact set in a metric 
Space. We will prove that the average error y, considered 
as a function of ay. is continuous on this compact 
set. Since any real valued continuous function defined 

on a compact set in a metric space has its infimum and 
Supremum on that set (37, 38) y has at least one minimum 
inside or on the sphere defined by equation (4.22). 
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y(C) is a continuous function on the compact subset defined 
by the above sphere. Hence y has an infimum value on this 
compact subset. This proves that in every sphere defined 


by equation (4.22) y has at least one minimum value. 
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In our case a = 0 and b tends to ~. The functions o,(t), 
i = 1,2,...N are, in general, oscillatory and cross the 
t-axis many times. The determinant D is zero if any col- 
umn of the determinant becomes zero. This happens when 
any function go; (t) has N zero crossings. Hence, in gen- 
eral, the Harr condition is not satisfied. 

However, we always approximate the semi infinite 
interval [0, ~] by a finite interval [0, Tol. Hence it is 
sufficient if the Harr condition is satisfied in this 
interval. The functions o; (t), dir= pli-2) s. 23N sane; Lincarly 
independent by definition. Hence if none of the functions 
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the interval (0, Tp] there exists a unique best approxima- 
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unique minimum lies outside this sphere. 


THE GRADIENT OF AVERAGE ERROR 

In the previous sections the existence and 
uniqueness properties of the minimum of y with a constraint 
on ISE are discussed. In order to compute the minimum 


point we can make use of the fact that at a local minimum 
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like o) (Ss). The minimization of average error with con- 


straints on ISE is discussed with examples in Chapter 6. 


BS. 
CHAPTER 5 
EVALUATION OF TRANSIENT RESPONSE 


INTRODUCTION 

This chapter develops a new numerical method of 
computing the transient response of a Laplace transform, 
expressed as a rational function of complex frequency. 
This method was primarily developed to compute the time 
domain representations of functions o(s), Kv" 72). 3 oN 
discussed in previous chapters. However this new tech- 
nique of computing the transient response is very general 
and computationally more efficient than some of the pre- 


vious methods (5). 
A BRIEF SURVEY OF PREVIOUS INVESTIGATIONS ON COMPUTATION 


OF TRANSIENT RESPONSE 


The evaluation of the inverse of a Laplace trans- 


form expressed as 


, a. # 0 (Sl) 


is commonly encountered in many branches of Engineering 
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and Science. The classical approach of evaluating the 
transient response is to express the rational function 
Y(s) as partial fractions and find the inverse of each 
factor. Generally this method is very complicated. More- 
over when two of the roots of the denominator are very 
close to each other, the evaluation of the corresponding 
residues of the partial fraction expansion leads to num- 
erical difficulties. Because of these problems, other 
numerical methods of evaluating the transient response 
have been developed. 

Corrington proves (40) that a rational opera- 
tional form of equation (5.1) leads to a linear n-term 
difference equation of the form given below 


y(t) = ? ne Be yy (t Sie et et (522) 


The coefficients FO 4 are real constants, independent of 
Ud 
t but dependent on At. These coefficients may be found 


from the following relation, 
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an easier way of determining the coefficients ae i= 
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Y(s) = 2 Stee (54) 
i=0 s 
The coefficients Cir ioe U2 en eare found) byAdividing 


the denominator into the numerator. If the Cc. coeffici- 


ents are known y(t) can be written as 


pig) Sige (5.5) 


Equations (5.2) and (5.5) are used to form the (n-1) simul- 


n-l 


taneous equations in (Fa) gel° 


Since Je is readily 
known from equation (5.3), these (n-1) equations can be 
solved to find aya Lbo= 1,2, <5. (n=1).. AFPurther |}computa— 
tions of y(t) are carried out by using equation (5.2). 
Aaron and Kaiser have pointed out some of the numerical 
difficulties one encounters while using this method (41). 
The division of two polynomials to find the Cc; coeffici- 
ents is difficult and leads to large errors (41). 

Another method of evaluating the transient 
response is reported by Liou (42). This method makes use 
of the state space approach. Liou considers the differ- 


ential equation which gives rise to the rational function 


W(s) of equation (5.1). ‘This™difterential~equation 1s 
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; ii 
where Dy (t) = oy (ty. The numerator coefficients of 


equations (5.1), Ge Daa may be expressed in terms of the 
initial values of y(t) and its first (n-1l) derivatives. 


Let Y(t) be a n-dimensional vector defined as 


T(t) =) (ei) Dye) whee Dey (E))) (5.7) 
Y(t) can be expressed as a state space equation 
Y(t) = A¥(t) (5.8) 
where T(t) = Sx¥ (t) and A is a nxn matrix. A is given as, 
CS me Ohana. mal 
0 0 5 Sa eas 0 0 
La i er? 
A = : , : : : (629,) 
1 0 i. ee Cie, |S 
a, “@n-1 ~*n-2 20 


The solution of the state space equation (5.8) is 


At 
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Y(t) = (07) (S010) 


where Y(07) is the initial condition vector. If we con- 


sider two points t = iT and t = (i + 1)T, we have 
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Equations (5.11) give, 
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Once the initial vector Y(o') is known, successive values 
of Y(iT) can be computed using equation (5.12). The ini- 


tial vector Y(0") is given by 


y(O" ) by 
eee y'(0") by - ary (0") 
(Om (5). 1'3)) 
n=2- + 
De aay (Om) : 1¥ (0 ) 
where yo (o* ) = dy (t) 
dt* ve 

t =O 


The transient response y(iT) is given by the first 
element of the vector Y(iT). 
One difficulty of the method is the necessity of 
evaluating the state transition matrix caer Liou suggests 
AT 


the power series evaluation of this matrix. e is expan- 


ded as follows. 
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(Sy. 4) 


where. Ac = [,wlcdentity matrix. 
‘ AT . : : : 
The’ matrix e is approximated by a matrix M, where 


Cayetissy), 


K 
e eM = 3 
i= 


0 


ge has developed a technique of choosing K 
such that the error involved in the approximation et as 
in equation (5.15) is within a specified limit. 

The computation of ere can be simplified by us- 
ing the relations between the last column and other col- 
ums of ae matrix. These relations are proved by 


Thomson (43) . 
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expressed in a recursive form as follows 
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modify the approach of Liou in order to avoid the compu- 
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where m._(t) is the element of Ae row and ae (last) 


shel 
column of Bee, From the definition of Sey in equation 


(5.7), we have 
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i=l 


In this method the CL values are computed first, using 
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m. (¢), 1 = 1,2;...,;n cam now be computed using equations 
(5.16) and y(t) is obtained from equation (5.26). This 
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NEW METHOD 
In the new method, y(t) is regarded as the out- 


put response of the system 


when an impulse input is applied. Let us consider a gen- 


eral system 


Y (s) 
Tien oe ee Cones) 


where Y(s) and U(s) are the Laplace transforms of the 
system output and input respectively. The system output 
may be evaluated by using the state space approach (45- 
46). The general differential equation of the system 


represented by the equation (5.27) is 
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substituting) y(t); (as giventby equation (5.29) into (5.28) 
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Equation (5.30) shows that equations (5.28) and (5.29) 


are equivalent. Hence 
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Css cle 3, 


The A matrix is the same as the one defined by equation 
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The equations (5.35) and (5.37) can be made use of to com- 
pute y(t) foneany ginpucru (tt). —Thessolution of "the state 


space equation (5.35) is given as (45-46) 
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y(t) = {Ce “Y(o) + Cfoe Bu(t)dt} + bou(t) Coa Oi) 


For the special case of the transient response of the sys- 
tem defined by equation (5.1) bo ==0'and= a Fey tan (eG the 
impulse at t = 0. The initial vector Y(o) can be consi- 


dered as that at t = 0. Hence 
Viole =" Viola or (5,41) 
This gives 
y(t) = cst eAlEM ag (yar 


Cea (peo) 


where 


(5543) 


The matrix A and vector B are given by the equations 


(5.9) and (5.35) repsectively. 


3 A : . = 
Let e : in equation (5.42) be approximated by 
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K terms. The approximate y(t) is then given by 


K-1 ,i,i 
y (Cb) =C. 82 ee. B 
c 2b 
1i=0 
DZ K-1, K-l 
7 Me FOE RS 
ae le) a 2! GF oOo om imeann tae o (5.44) 
: = 


Since factors like = are scalars y(t) can be written as 


2 K=1 
Y(t) = Cele ae aet a5 — cage suet —_ SEs (5245) 


The matrix inside the bracket in equation (5.45) can be 
expressed as the product of a partitioned matrix and 


another vector. Thus we get 


e 
K-1 
, enna te Meer een t 2 t | 


Since the vector B has zeros for the first (n-l) elements 


and unity for the last element we have 
[rs iapjatsi-.----ia-2p| = E ety eee A |- ® 
(5.47) 


where A; is the last column of A” Matrix, A’ = tee the he encd cy, 
Matrix and © is a nxK matrix. The ios column of the 6 
mMAtEix is theelast colummeor the siatriseA ew) Hence there 


exists a recursive relation between any two successive 


columns of the @ matrix. This relation is proved as fol- 
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lows. 


AY Gena (5.48) 


i 


Let the last column of A” ~ be given in terms of the ele- 


ments of the $@ matrix. That is 


c 


Sale Meas Lee. cee 


h th 


where da% is the element of 4° row and k column of 


matrix. 
From equation (5.48) we get 


A. 


A A. 


a i-l 
That is 
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(5250) 
The evaluation of equation (5.50) gives the recursive 


relation as 


awol 


en? 36d 


Ss z5iw 


(we8Mm 


err 


— 
ig Moly 


- 
+ 
fr 
4 
~ 


98. 


Teak oom 

Ere ae 
|= | | (5.51) 
Pe Oe ieee | 

acu : pepe neitl yi 


This recursive relation can also be expressed as 


4 it] = Pasa il SS Lo yee ne) (S52) 


n 
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They relation (5252)" is valad for all valves of i% “tence 
knowing the first column of 6, which is the last column 
or Ace tiateis, the Last column of the Identity matrix, 


the transient response can be computed from 


y(t) = Coy (5.53) 


2 K-1 Jt 
where v = E cr Sr aon f 


(5.54) 
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Using equations (5.53) and (5.54) y(t) can be expressed 


as 


K 2 
y(t) = rh.t (5.55) 


pee oe (5.56) 


> 


where h; = . (Db; 
j=1 

The evaluation of each column of $ matrix and the corres- 
ponding h, coefficient may be done simultaneously. Once 
the h; coefficients are computed and stored the transient 
response y(t) can be evaluated for any t. 

The error involved in approximating y(t) in the 
above manner is now studied. We make use of the approach 
of Liou (42) to find an upperbound of the error. The 


exact expression for y(t) is 


y(t) = ceMtp 

a ob pal 
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We’ are approximating y(t) by the first term of the right 


hand side of equation (5.57) 


y(t) =C B (5258) 
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The error ‘e(t)® 1s"given by 


co ieee 
a(t) =6C) ee eoeeE (5.59) 
: ue 
1=K 
= CRB 
ies ries 
where R= 2 Ba (5.60) 
i=K a 
Retonacnixnamat lin, 
Mg Cmax 1S an upper bound of each element ee of R so 
that 
Wee IM k4 hs (5262) 
iy max 
we obtain from equation (5.59) 
n 
pone) Soe ole: (5,62) 


An upperbound rn can be found by defining a suitable 


ax 


norm for the matrix A. Liou suggests the norm as 


m n 
Mall = 2 jo,,| (5.63) 
1,j=L 
where oi 5 is the ‘capes element of A matrix (42). Later 
he has modified this definition of ||A|| (46) as, 
|) (5.64) 
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It is found that in most cases K, the number of 
terms estimated by using the definition of ||A|| in equa- 
tions (5.63) and (5.64) is much larger than what actually 


is required. Hence we make use of a modified definition 


of the norm of A as 


[Jal | = max (2 [one |) (5.65) 
a (ME 


Fe 
j J 
According to this definition, the norm of A matrix of 
equation (5.9) is 

PEM eee ese) SNR [ro abe a) (52°66) 

i=l,n-1l 

Tee la, | >> 1, then the norm of A is equal or approximately 
equal to the largest magnitude of ass Paw 2c pti, We 
now prove that the definition of equation (5.65) satis- 


fies all the conditions of a norm. These conditions are 


1) A # 0 implies ||A|| > 0 
2) |[aal| = |a] |{A]|, o is a scalar 
2) ee 
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se See a eo 


where A, and ae are any two (nxn) matrices. 
By definition (2S) nee are not zeros. Hence con- 


dition 1 is always satisfied. 
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BOonsCondi tions jp alet 
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(57,08) 
This proves Condition 2. 


Condition 3 is proved as follows. Let O, and Bis be the 


J 
Gugee element of the matrices A, and A, respectively. 
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A fourth condition necessary in finding an upperbound of 
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EQUaCION (5.1/3) ells aceuc]toOreanyece Lunn Or git and hence 
wergec 
k+1 k 
HAT} s TAT] LAT 
Sola zeae (5.74) 


Applying the above properties of the norm to equation 
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Thus we get an upperbound on e(t) as 
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By making use of equation (5.79) it is possible to choose 
K such that the error is within allowable limits. 

This method of computing the transient response 
has many advantages. It is not necessary to compute as 
Matrix and the initial vector as required in Liou's method 
(42). The method suggested by Valand involves the manipu- 
lation of n series while the new method involves only one 
series for all values of n. The distribution of the 
eigenvalues of the system does not affect the method. For 
example, even if two poles are identical or very close 
to each other the method is capable of giving the response 
without any further modification. Since y(t) is expressed 
aS a power series in t the derivatives at any value of t 
are easily obtainable. The round off errors of computa- 
tion do not propagate with increasing values of t. This is 
because y(t) is evaluated at each t by the power series 
and hence is exact. 

However, when ||A|| is very large and it is 
required to compute y(t) at a large value of t some dif- 
ficulties are encountered. Under such circumstances it 
is necessary to consider a large number of terms. One 


method of solving this problem is to split Y(s) into 
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smaller factors. The denominator of these factors of 
Y(s) can be formed from the poles of the system. The 
numerator coefficients of the factors of Y(s) may be 
found by solving the n simultaneous equations obtained by 
substituting n values of s. But a better way of doing 
this is to form n simultaneous equations by using the 
initial value of y(t) and its first (n - 1) derivatives. 
These values are easily obtained by finding the first n 
coefficients of y(t) as discussed in these chapters. 
These coefficients are equated to the corresponding coef- 
ficients of the factors to obtain the n simultaneous 
linear equations. The solution of these equations deter- 
mine the above factors. As ||A|| of each of these fac- 
tors has smaller value, the number of terms K necessary 
in each case is also smaller. The transient response 
y(t) is then obtained as the sum of the transient responses 


of the individual factors of Y(s). 
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CHAPTER 6 
COMPUTER ALGORITHMS AND NUMERICAL EXAMPLES 


INTRODUCTION 

In this chapter we first discuss the computa- 
tional techniques employed to find the best pole posi- 
tions by minimizing the ISE. Instead of employing any of 
the iterative techniques discussed in Chapter 3, we rely 
on efficient numerical methods of minimizing functions of 
several variables without calculating derivatives (47, 
48). This avoids the complicated filtering operations 
required in all the exact methods of minimizing ISE 
reviewed in Chapter 3. The minimization of average error, 
with constraints on ISE, is done by using the penalty 
function approach (51). Finally the method is applied 


to a specific example which is an ideal low pass filter. 


COMPUTATION OF ISE 
The concept of the complimentary filter is used 
to’ compute the ISE for any set Of pole positions. The 


expression for the ISE, aS given in Chapter 2, is 
° i= 2 
ISH = I = fj) a(t) |odt (6,1) 


where a(t) is the output of the complimentary filter when 
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Vit)iaethne reversed f(t), is appldedvas an input «This 
relation was discussed in detail in Chapter 2. The com- 


plimentary fidter Gi(s) .1s. given. by 


The scheme for evaluating the ISE is shown in Fig. 2.1. 
The given function f(t) is time reversed and applied to 
the filter G(s). Since any physical signal vanishes for 
some t > Tos the reversed signal v(t) may be considered 

as starting at t = 0 and extending to t =T,. The expres- 


sion for the ISE—-then becomes 


lat (6.3) 


Ise = 1 = se°|a(t) 
The output a(t) of the complimentary filter is computed 
By wmsing equation (5.36) r8*Thes equatrzon= givess the general 
expression for the output of any filter’G(s) which has 
the form given by equation (6.2). Thus the expression 


for+atty Zs 


Equation (6.4) involves the evaluation of the state vari- 
These variables are defined by equations 


(S70. = 5.34). The solution of equation (5.34) gives) che 
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tial values ‘of Y;, are zeros. 
is computed by evaluating the 


The integration uses the Five 


A subroutine ERROR is written 


routine accepts a set of coefficients (a; ) 
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This is accomplished by the Runge-Kutta 


a standard subroutine in 


Package “Libracy /?*rhe ini 


Once a(t) is known the ISE 


integral of equation (6%3):. 


point Ouadratureyrormula. 
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i=l! the reversed 


Signal v(t) and computes the corresponding ISE of least 


square representation. 
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MINIMIZATION OF ISH 


This program is given in Appendix 


The general approach employed in minimizing the 


ISE is to consider the ISE as 


a function of the parameters 


Giachne complimentary filter and to find) these paramecers 


such 


that “the ISE «2s minimized; 


This has the advantage 


Of reddciig “the original 2N variables of Algrain and 


Williams equations to N variables of the complimentary 


filter. Another advantage of 
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are considered as the variables, 


this "approach “is “that, "ar 


the operations 


involving complex numbers can be completely avoided. 


A direct method 


one of the gradient techniques. 
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pole positions over a wide range near the optimum point. 
This makes the gradient techniques very inefficient in 
Minimizing the ISE. Since we want to avoid the compli- 
cated filtering operations necessary for the various 
iterative schemes, a new method of minimization was tried 
and found to be very useful. This method does not require 
the computation of the gradients. The basic approach 
employed is due to Powell (47). Powell uses the method 

of conjugate directions to minimize functions of several 
variables. Zangwill has pointed out (48) certain draw- 
backs in Powell's method and has suggested a new algorithm. 
This algorithm incorporates all the basic features of 
Powell's method but does not suffer from the difficulties 
pointed out by Zangwill. We have made use of the Zangwill 
algorithm with some modifications, to minimize ISE. The 


basic theory of these algorithms is now discussed. 


FUNCTION MINIMIZATION WITHOUT CALCULATING DERIVATIVES 
(1) Powell's Method 

The minimization of a function of several vari- 
ables without calculating derivatives has attracted the 
attention of several investigators (47-51). Among these 
the method of Powell (47) is considered the best as this 
will minimize a quadratic function of n Variables in n 


iterations. Powell considers the minimization of a quad- 
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Let 


Ox) = XxX PX eB XK ec (655) 


where P is a nxn matrix, B- is a n- vector and c a scalar 
constant. Let (Ea) ey be n conjugate directions such 


that 
2apey = hye tee (6.6) 
Le) : 


Let Xo be an initial point. Any point X may be expressed 


in terms of these conjugate directions. 


Let 
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Simplifying by using equation (6.6), we obtain 
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Q(X) = Q(X,) + 2 oi OPE + E;PX,) 


poe (6.4) 
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Examining equation (6.7) it is found that minimization of 
Q(X) by searching in one of the conjugate directions is 
independent of the other conjugate directions. Hence 
minimum of Q(X) may be found by searching along each of 
Ehesconjugatesdi rections only oncey The problemiot mania 
mizing a quadratic function is thus reduced to the prob- 
lem of finding n conjugate directions. Powell suggests 

a method of determining these conjugate directions. Let 


Xo and X, be two points such that the quadratic Q(X) is 


L 
minimum at Xo and X) along a direction jn. This means 
99 (x + an) = 0 ata= 0 
eRe) 0 
0 = = 
and Fae (Xy “Cir ct)) ee eel tae 0 (678) 
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Substituting ® = 0 and using equation (6.8) we obtain 
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t 2 
n P(X) = Xo) = 0 (Ong) 


By equation (6.11) the direction (X) - X54) is P conjugate 
to n. Based on this theory Powell suggests an algorithm 
to minimize Q(X). 
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jizz 2%¢ chosen as the co-ordinate dir- 


ections even though they are not P conjugate directions. 
At the end of n searches along these directions a point 


xX) is obtained. A new direction Ent] is chosen such that 
=X, - X (on 22) 


The function is minimized along ntl COnyteld asnew ~start— 


ENGADOGINe eX The directions are rearranged according to 
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the rule 


ay = Eas Sn (65.153)) 


The iteration is repeated until the function values does 
not change. 

The function Q(X) will be minimized inn itera- 
tions is proved as follows. At the end of n linear 


searches of each iteration a new direction is chosen by 
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using the relation (6.12). The function Q(X) is minimum 
at xX) and Xo along the direction BAe Hence the new direc- 
tion chosen at the end of any iteration is P conjugate to 
En’ Thus at the end of n iterations the function Q(X) 
has been minimized along n mutually conjugate directions 
and so the minimum of Q(X) has been obtained. 

One drawback of the method is that the new direc- 
tions chosen need not be linearly independent. The n 
direction generated need not span the n dimensional space, 
in which case the method will not give the minimum of 
Q(X). Zangwill has shown with a counter example that 
Powell's method does not converge in the case of a partic- 
ular function which is strictly convex, quadratic and has 
a unique minimum (48). This led him to suggest a new 
algorithm. This new algorithm with some modifications, 


is now discussed. 


(2) Zangwill's Method 


This method incorporates the quadratic conver- 
gence properties of Powell's method and guarantees that 
in the case of a quadratic function the new directions 
chosen are linearly independent. This is achieved by 
ensuring that any new direction chosen can never be zero. 
The algorithm, as modified to suit our problem, is as 


follows. 
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Let Pan be the n normalized co-ordinate 
aLrrections, and oven be the yn normalized conjugate direc-— 
tions used in the rae iteration. initial conjugate direc-— 


tions are chosen as the co-ordinate directions. Xo is 
n L 
the initial point and ||[x|| = ( 2 cone is the usual 
i=l 
Eucledian norm. It is assumed that P matrix is positive 


definite. (xt) denotes a point obtained by minimizing 


along mage conjugate direction in Ree iteration and (ay 


th 


denotes the starting point of K iteration, 


INITIALIZATION 
; ) i: 
Let us define x ad and Xp as 
x? = ee = X iniveral point 
nt+1 0 Oy ; 
ote ‘ (6.14) 
aL ae 4’ pe oe ° 
K = 1 
and 4 = 1, where K and 3} are indeces. 


(a) Minimization along conjugate directions. 
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The new conjugate directions for the next iteration are 


chosen as 


2 le Ar oly (6.128) 


(b) Minimization along co-ordinate directions. 


Choose a to minimize yee ‘oar. ). ff a =f 


i 3 
and j =n, repeat (b) with j = 1. If a = 0 and j #n, 
Feneat (5) with y= j +i. j&-Lf {(b)° is repeated n times 


in succession sare LS a Minimum point. If a 4°90, choose 
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Parts ta) and (b) of iteration Kare repeated in succes— 
sion until the minimum point is obtained. 


In the above algorithm the function is first 
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minimized along the n conjugate directions and then along 
one of the co-ordinate directions. If the function does 
not change along this co-ordinate direction the next co- 
ordinate direction is tried. If all the n co-ordinate 
directions are tried in succession and there is no change 
in Q(X), Grad(Q) at this point is zero and hence a local 
Minimum is obtained. The major improvement of this method 
over Powell's method is in the choice of the new direc- 
Lon an in step (a) of iteration K as given by equation 
(6.16). This ensures that bea is always different from 
zero, except for iteration l. If a = [0], it means 


that the initial point is a minimum point and the search 


can be immediately terminated. For K > 1 we note that 


K K he 
OUR OK 1 gale a) (6.20) 
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and hence a x X The point a is a minimum along ES 


and the point X is also a minimum along the same direc- 
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tion. Therefore aan is P conjugate to aw Since Enel 


and ES are non zero vectors and the matrix P is positive 
definite, it follows that these directions are also lin- 
early independent. Hence the quadratic function Q(X) can 
be minimized in n iterations. 

It can also be proved that if the function Q(X) 


is. strictly convex, is continuously differentiable, for 
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all X and has a greatest lower bound, the method converges 
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The method guarantees that 
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Since Q(X) has a greatest lower bound, the relation (6.21) 


gives 
lim Q(x**+) 
K>o 
Let us define 
lim Q(xt1) 
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and Lim O¢ xX.) 
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where eee is some conjugate direction and o(x* + an) is 
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Since Q(X) is continuously differenti- 
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Q(R) = Q(X + an) < Q(R + Bn) (6.24) 


where 8 is any scalar. 


Since Q(X) is strictly convex there should be at 


least one point between % and Rais along ny such that 
+1 


at this point x, O(x) 1s Jess than Q(z) or Q(z ). As 
this is not possible, 
at tnieg® (6.25) 
The point eS is arrived at, after competing step (b) 
of iteration K of the algorithm. Hence 
+1 00 
Q(z PESPO UR? + Br5), 3 eae Sa Sa (6526) 


Equation (6.26) ensures that the final point x is an opti- 
mum point. 

The original method of Zangwill requires the 
availability of n normalized directions different from co- 
ordinate directions to start with (48). The method as 
given here does not require this. The algorithm has been 
suitably modified so that the co-ordinate directions alone 
are sufficient to start with. 

A digital computer program ZAGMIN is written in 


FORTRAN IV which will minimize any function of several 
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variables using the above method. It is found that the 
criteria for terminating the minimization as suggested in 
the algorithm is too strong. Hence the program is also 
Made to terminate when the change in function value in 
successive iterations is less than a preassigned value. 
The initial step of the linear search is to be specified. 
This step is doubled or halved at each search point until 
the minimum is passed. Three points at equal intervals 
between which the minimum lies are found and the minimum 
point is computed by quadratic interpolation (47, 51). 
If the function value at this point is not less than that 
at the middle point, the computed point is ignored and the 
middle point is taken as the minimum. 

The ISE, I, is considered as a function of the 


variables (a and is minimized for different functions 


i) i=l 
using the above program. It is found that the new method 
is more efficient than the original method of Powell in 
Minimizing the ISE. The values of the parameters (ee) es 
at the minimum ISE give the best pole positions. These 
poles are used throughout the subsequent analysis and 
design. 

Once the poles are determined, the orthonormal 
functions o(s), kK = 1/2, acc can be found... The order of 


the poles is not important in forming these orthonormal 


functions. In the examples considered the real poles are 
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taken first, followed by the complex pole pairs. A digi- 
tal computer program for computing the numerator and 
denominator coefficients of the orthonormal functions 
(Ss), k = 1,2,...N, is given in Appendix I. The inverse 
of these functions d,(t), k = 1,2,...N are computed using 
the methods discussed in Chapter 5. The subroutine 

SERIES takes in the numerator and denominator coefficients 
Ofmany  Gatronal ruiceion o) (s) and computes the coeffici- 
ents of the series expansion of , (t) in powers of t. The 
maximum value of t at which p; (t) is required is also given 
aS an input to SERIES. If this maximum value ty > 1, the 
series expansion is made with respect to a new variable tT 


such that 
T = t/t (Gu2) 


The total number of terms necessary can be calculated by 
using the relation (5.79). Another method of determining 
the number of terms is to make use of the condition that 
the absolute value of te term of the series expansion 
tends to zero as n tends to infinity. As all computations 
are carried out in double precision arithmetic the compu- 
tation of the coefficients is terminated when the terms 
16 


become less than 10 °°. The relation (5.79) is used as a 


general guide to find the total number of terms required. 
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The function subprogram VAL uses these coefficients and 
maximum value of t to compute the inverse at any point t. 
The subroutine SERIES and function subprogram VAL are 


given in Appendix I. 


MINIMIZATION OF AVERAGE ERROR 
It has been proved in Chapter 4 that the fre- 
quency response of H(s) can be improved by minimizing the 


average error. Hence the approximation 


N 
h(t) = 2 C.9¢, (t) (6.28) 


is to be so chosen such that the average error y is mini- 
mized subject to the condition that ISE is less than a 

preassigned value yw. This can be achieved if the minimi- 
N 


Zation OL ¥ 12S done such that (C5); are always chosen 


=] 


to satisfy equation (4922). This equation is 


N % av . 
a ape) ae eg eart (6.29) 


where 
Cae fof(t)9o, (t)dt 
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This constraint on the choice of oe can be incorporated 


into the expression of average error by the penalty func- 
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tion approach (51). The expression for the average error 


is modified as follows 
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Yn = folf(t) - 2 
(6.30) 


where w is a weight factor and p is a constant defined as 


N % 
Ort 0 when 02 (Come —6 C28) =a 


(6.31) 
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The value of the weight factor w may be changed as desired 


(51). The choice of p as in equation (6.31) ensures that 
when (ae is inside or on the sphere defined by equa- 
ieLOne (0.29) > tes and y are the same and whenever the search 


point goes outside this sphere, van is made greater than y. 
The actual value of ies is controlled by a proper choice of 
w. it has been proved in Chapter 5, that y has a minimum 
value in the sphere defined by equation (6.29). Hence as 
the minimization of i is continued the points of search 
approach this sphere and finally converge to a point within 
Or On this sphere. 

Any of the gradient techniques may be used to 
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* N puso? Que 
Grad Tew, = Grad (y) + wpGrad| © (Cc. - C.)*~ - 6 | 
: an i 
i=l 
(6.32) 
th ; 
The k component of Grad Ope? is given as 
2Yn By v ¢ N v2 2 
ma = = + 4wo(C, - C_)¥ XY (C. - C.)” - 6 } (6.033) 
dC) dC, k ke ie i 
i=l 


The term am is computed as discussed in Chapter 4. This 
k 
Ss, 


is given a 


ro = fa°Sgn{e(t) }{-$, (t) }at (6.34) 
N 
where ea(e=— £(t) — 2 C. o>, (t) 
ail 


The interval [0, To] is divided into small intervals and 
the integral in each interval is evaluated using Gauss 
quadrature formula. The error involved in evaluation of 
the gradient due to the discontinuity of Sgn(e(t)) can be 
made negligible by choosing the subinterval to be very 
small. This avoids the determination of the zero cros- 
sings of e(t) for each set of coefficients conan Simi- 
larly the time of computation of average error and gradi- 
ent is considerably saved by computing and storing the 
values of the functions o,(t), 1431, 2,0.0,N at. thOSe. pointes 
as required by Gauss quadrature formula. This need be 


done only once at the beginning of the minimization. The 
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subroutine MDQG makes use of these values and computes 
yN 
i‘i=1° 


numerical examples the actual minimization of average 


ie and Grad ye? for any point (C In the following 
error was carried out by using Fletcher-Powell algorithm 
(52). 

The subroutine MDQG and a program for minimiz- 
ing the average error are given in Appendix 1. The mini- 
mization of average error makes use of the standard 
Pleecher=Poweli “algorithm available in the IBM Scitentiiae 


Subroutine Package Library. 


NUMERICAL EXAMPLES 

The method developed in this dissertation was 
used to approximate the ideal low pass filter. The rea- 
son for choosing this example is that it can be made use 
of to check the comparative merits and demerits of the 
new method. 

The frequency response of an ideal low pass is 


defined as 
F(jw) = e / 0, Jol <1 
= (0) | w| Su (6.539) 


The corresponding time function f(t) is given as 
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f(t) is obtained by delaying the sint/t function by to: 
These functions are shown in Fig. 6.1. 
Theoretically F(jw) is not realizable. This is 

because f(t) exists for negative values of t. The func- 

tion f(t) is not absolutely integrable. Hence f(t) does 
not satisfy the assumptions made in Chapter 4. But even 


though F(jw) is not physically realizable, it is possible 


to obtain a satisfactory approximation of F(jw) by suit- 


ably truncating f(t). Two different versions of f(t) are 
used to approximate F(jw). They are 
Peles Lota a) 
Dee es) = a gates , cCepoy, 3m.) 
= 0 PRCe Oy. oT (Gna) 
sin (tyes 0,) 
2) eo) = iT, MESA at te 10, 47) 
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f(t) as defined by equations (6.37) and (6.38) are physi- 
cally realizable. These functions are shown in Fig. 6.2a 


and Fig. 6.2b respectively. Each of these functions is 


h h 


approximated by 5 order and ai order filters. These 


approximations are now given. 
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Figure 6.2. Truncated Impulse responses of ideal low- 


pass filter. 
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The following notations are used in the results. 
~ 


H(s) - Least square filter. 
H(s) - Filter obtained by minimizing average error y 
with constraints on ISE. 
ING N AW) 
(bs) soy - The numerator parameters of H(s). 
Crea - The numerator parameters of H(s) 
NG 
(an) =; - The denominator parameters of H(s) and H(s). 
NG 
(C;) - The coefficients of orthonormal functions at 
Minamum ILSh. 
cepa - The coefficients of orthonormal functions 
at minimum average error with constraints 
on ISE. 
N AV) 
(S. 52) - The poles of H(s) and H(s). 


av) av) 
h(t) and n(t) are the time functions of H(s) and 


H(s) respectively. 


@ (t) and e(t) are the corresponding error functions. 


The orthonormal functions for each case are 


given in Appendix II. 
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a) Order of filter, N = 5 
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The frequency responses of H (ju) and H(jw) are given in 
Fig. 6.3 and Fig. 6.4 respectively. The least square 
approximation h(t) and the error function @ (t) are given 
in Pig. 6.5.) Flo. 6.6 gives the function (ct) and the 
corresponding error function e(t) for minimum y with con- 


straints on ISE. 
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Fig. 6.3 Magnitude and Phase responses of H(jw), N=5, Delay=n, Minimum ISE=0.00021. 
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Fig. 6.4 Magnitude and Phase responses of H(jw), N=5, Delay=n, Minimum y=0.06163, ISE<0.0005. 
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Figure 6.5. The least square approximation A(t) and error e(t). N=5, Delay=n, Minimum ISE=0.00021 
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Figure 6.6. The approximation h(t) and error e(t) for minimum average error. N=5, Delay=n, Minimum y=0.06163, ISE<0.0005 
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b) Order of filter, N = 8 
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The frequency responses of H (jw) and H(jw) are given in 
Fig, 6./ and=Figq., 646 respectively. “The least square 
approximation h(t) and the error function € (t) are given 
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Figure 6.8 Magnitude and Phase responses of H(jw), N=8, Delay=7, Minimum y=0.02058, ISE<0.0005. 
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Figure 6.9. The least square approximation R(t) and error e(t). N=8, Delay=n7, Minimum ISE=0.000057 
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Figure 6.10. The approximation h(t) and error e(t) for minimum average error. N=8, Delay=n, Minimum y=0.02058, ISE<0.0005 
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the function ir(t) is given in Fig. 6. 2b. 
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The frequency responses of H (jw) and H(jw) are given in 

Fig. 6.11 and Fig. 6.12 respectively. The least square 

approximation. h(t) and the error function @ (t) are given 
an Fag. 6.13.) Fig. 6.14 gives, the function.) Ge) and ene 
corresponding error function e(t) for minimum y, with 


COnNStraints on ISE. 
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The frequency responses of H (jw) and H(jW) are given in 


Fig. 


Oielo and igi. 


6.16 respectively. 


The least square 


approximation h(t) and the error function 8 (t) are given 


THe ED, Ocul. 


corresponding 


Fig. 


6.18 gives the function h(t) 


and the 


error function e(t) for minimum y, with 


constraints on ISE. 
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SUMMARY AND CONCLUSIONS OF NUMERICAL EXAMPLES 

The numerical examples show that it is possible 
to improve the frequency response by minimizing the aver- 
age error. In examples 2a and 2b there is considerable 
improvement in the frequency response while in examples 
la and lb this improvement is not significant. In example 
1 it was possible to reduce the ISE to a very low value 
and hence the unconstraint minimum of average error was 
very close to the centre of the sphere of equation (6.29). 
Therefore the minimization of average error gave a set of 
coefficients very close to those at minimum ISE. But in 
example 2, it was not possible to reduce the ISE to such 
low values as in the first example. The coefficients 
giving the unconstrained minimum of average error, were 
not close to the centre of the sphere of equation (6.29). 
In this case it was possible to improve the frequency 
response by minimizing the average error while keeping 
the ISE less than some suitable value. It was observed 
that the magnitude characteristics of these improved 
responses had a constant attenuation at the low frequency 
(w < 1). This was due to the relaxing of the condition 
on ISE. This could be overcome by any one of the fol- 


lowing three methods. 
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(1) By reducing the radius of the sphere giving the con- 


straint on ISH; that is, by reducing ISE. 
(2) By multiplying H(s) with a suitable constant. 


(3) By keeping the constant term of the numerator of 


H(s) and H(s) to be the same. 


In the examples 2a and 2b the frequency responses of H(s) 
were adjusted by method (3). It was found that this gave 
the desired response at low frequencies without affecting 
the improvement obtained at higher frequencies. 

The frequency responses of these filters at higher 
frequencies are inferior to those of the corresponding 
Butterworth filters. But, these filters have better 
frequency characteristics in the transition band. For 
example in case 2b (N = 8, Delay = 27, y minimum) the 
maximum ripple in the pass band has a value of -43db. 

THLS Magnitude 1s first attained atiw = 1554. “The Butcer- 


worth response of ou order filter at w = 1.54 is only 
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CHAPTER 7 
CONCLUSIONS AND RECOMMENDATIONS FOR FUTURE WORK 


The problem of approximating an arbitrary fre- 
quency response by a realizable filter has been studied 
in this research. The general approach employed was to 
convert the frequency domain approximation into an equi- 
valent time domain approximation. This was done by mak- 
ing use of two important relations between the frequency 
domain and the time domain errors. These are the equiva- 
lence of least square criterion in the frequency domain 
and in time domain and the relation between the least 
upper bound of the error in frequency domain and the aver- 
age error in the time domain. A new method of minimizing 
the upper bound of the magnitude of the frequency domain 
error with constraints on integral squared error has been 
developed. This method made use of a set of orthonormal 
functions of exponentials. 

One important aspect of the problem was the 
determination of the poles of the filter which would yield 
the least integral squared error. The concept of compli- 
mentary filter was used to evaluate the ISE. This has 
the advantage that the ISE corresponding to any set of 


distinct poles could be easily computed using a simple 
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filtering operation. This did not involve computations 
using complex numbers. An algorithm was developed and 
implemented to compute the ISE. The minimization of ISE 
was carried out by using Zangwill algorithm. This algo- 
rithm was an improvement on Powell's method of minimizing 
a function of several variables without calculating deri- 
vatives. The ISE is very insensitive to variations in 
pole positions over a wide range around the optimum point 
and hence the gradient methods of minimization was not 
successful in minimizing ISE. This led to the use of 
Zangwill algorithm. This algorithm with some modifica- 
tions was implemented. The mathematical basis of this 
method was also reviewed. 

A new method of choosing the zeros of the fil- 
ter so as to minimize the deviation in the frequency do- 
Main was developed. It was proved that the average error 
in time domain is related to the least upper bound of the 
magnitude of the deviation in the frequency domain. This 
gave a convenient way of minimizing the deviation in the 
frequency domain. However, the minimization of the aver- 
age error does not mean the minimization of ISE and vice 
versa, It was found that minimization of one of these 
errors alone need not give the desired result. This led 


to the development of a method of minimizing the average 
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error subject to constraints on ISE. This was carried 
out by using a set of orthonormal functions of exponen- 
tials. The mathematical basis of this scheme of minimiz- 
ation was established. It was proved that for every ISE 
less than or equal to a preassigned value there exists 

a set of coefficients of the orthonormal expansion such 
that the average error in time domain at these coeffici- 
ents is a minimum. The uniqueness properties of this 
approximation were also discussed. The actual minimiza- 
tion was carried out by using the Fletcher Powell algo- 
GithMaeethe constraints onsthe ,ISkeweresimplementedsby 
using the penalty function approach. 

In order to implement the above scheme of mini- 
mization it was found necessary to compute the inverse of 
a Laplace transform expressed as a rational function of 
complex frequency. A new method of numerically computing 
this inverse was developed. In this method the transient 
response was expressed as a serieS expansion in powers 
of the independent variable t. The coefficients of the 
series expansion could be easily computed by using the 
recursive relations established for this purpose. The 
new method was found to be more efficient than some of 
the already existing methods of evaluating the transient 


response (5). This method makes the approximation of an 
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arbitrary frequency response using the time domain approach, 
more practical. 

The results obtained by applying this technique 
to the specific example of approximating the ideal low 
pass were very encouraging. The frequency characteristics 
of the filters obtained by this method were better than 
those of the least square filters. This improvement in 
the frequency response is not very marked in the first 
example considered (Delay = 1) while in the second example 
(Delay = 27) the improvement was very significant. This 
was because of the fact that in the first case it was 
possible to minimize the ISE to a very low value and 
hence the two points of minimum ISE and minimum average 
error were close to each other while in the second case 


this was not true. 


RECOMMENDATIONS FOR FUTURE WORK 

In the example considered, the time domain 
function of the desired frequency response was analyti- 
cally known. In general, when approximating an arbitrary 
frequency response the corresponding time domain function 
should first be computed. This can be done by using the 
Fast Fourier Transform algorithm. It is also possible to 


check beforehand by using the inverse FFT whether. the 
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truncated time function can give a sufficiently close 
approximation of the desired frequency response. Once 
the time function is known it should be reversed in time 
before applying as input to the complimentary filter for 
evaluating the ISE, 

The extension of the method to the design of 
Recursive Digital filters for arbitrary frequency response 
Offers an interesting field for future work in this direc- 
tion. A set of discrete orthonormal functions have al- 
ready been developed by Young and Huggins (54) and Broome 
(55). The least square error of approximation can be 
evaluated by using the concept complimentary filter in 
Z-domain (56). All the theories developed in this dis- 
sertation can be extended to the discrete system. The 
new method of inverting the Laplace transform can be modi- 
fied to develop a corresponding numerical method invert- 
ing Z-transform. The method of computing transient 
response can also be made use of to find the equivalent 
Z-transform expression of a Laplace transform (57). 

Another field of future work is to use the 
optimum poles to construct a set of functions satisfying 
the Harr condition and to carry out, the minimization Of 
the average error with respect to the coefficients of 


these functions. This will improve the quality of the 
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approximation, 

This research has shown that it is possible to 
get a better frequency response by relaxing the condition 
on ISE. This suggests the possibility of directly mini- 
mizing the Chebyshev error in the frequency domain without 
using the relation between the upperbound in the frequency 
domain and the average error in time domain. This involves 
the evaluation of an analytical expression for Chebyshev 


error in frequency domain. 
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APPENDIX I 


SUBROUTINE ERROR (NDIM,AC,VALUE,Y,DERY, AUX, RFX) 
EXTERNAL FCT,OUTP, RFX 
SUBROUTINE ERROR (NDIM,AC, VALUE, Y,DERY ,AUX, RFX) 
THES SUBROUTINE COMPUTES THE INTEGRAL SQUARED ERROR 
(ISE) OF EXPONENTIAL REPRESENTATION BY THE METHOD OF 
COMPLIMENTARY FILTER.IT USES THE RKGS SUBROUTINE IN 
SSPLIB OF IBM.THE ARGUMENT LIST OF RKGS IS CHANGED BY 
ADDING A NEW PARAMETER AC OF N VECTOR.AC IS THE VECTOR 
CONTAINING THE COEFFICIENTS OF THE COMPLIMENTARY FIL 
TER.A(1) IS THE COEFFICIENT OF S**(N-1) OF DENOMINATOR 
OF COMPLIMENTARY FILTER. 
THE COMMON STATEMENT IS USED WITH SUBROUTINES ERROR, 
ZAGMIN, -ECT,AND OUTP 
NDIM # OF THE COEFFICIENTS 
AC-NDIM VECTOR CONTAINING THE COEFFICIENTS OF DR OF 
COMPLIMENTARY FILTER 
VALUE-THIS GIVES THE ISE ON RETURN 
Y-VECTOR OF DIMENSION NDIM. CONTAINING THE STATE 
VARIABLES.INPUT VALUES OF THIS VECTOR ARE ZEROS. 
DERY-VECTOR OF DIMENSION NDIM . INPUT VALUES ARE 
SPECIFIED BY RKGS. 
AUX-AS SPECIFIED BY RKGS 
RFX- A FUNCTION SUBPROGRAMME GIVING THE TIME REVERSED 
SIGNAL 
DELT-INTERVAL OF INTEGRATION.ERR- ISE. V-VECTOR OF 
DIMENSION 5,A WORKSPACE. 
BV,S,KOU ARE WORKSPACES. 
TL-DURATION OF SIGNAL 
REAL PRMT(6),Y(10) ,DERY (10) ,AUX(8,10) ,AC(10) ,V(5) 
COMMON DELT,ERR,BV,V,S,TL, KOU 
ERR=0.0 
KOU=0 
BV=0. 
S=4,*DELT 
PRMT (1)=0. 
PRMT (2)=TL 
PRMT (3) =DELT 
PRMT (4)=1.E-4 
PRMT (6) =DELT 
X=NDIM 
0: on ny p.4 
DO 7 I=1,NDIM 
¥(1)=0. 
7 DERY (I) =xx 
CALL RKGS (PRMT ,Y,DERY,NDIM,IHLF,FCT,OUTP,AUX,AC) 
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VALUE=ERR 
RETURN 
END 


SUBROUTINE FCT (X,Y,DERY,NDIM,AC, RFX) 


die 27 


C THIS SUBROUTINE COMPUTES THE DERIVATIVES OF STATE 
C VARIABLES OF COMPLIMENTARY FILTER AT X. 
C ALL PARAMETERS ARE AS GIVEN IN ERROR. 
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& 
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ALL PARAMETERS ARE AS GIVEN BY ERROR AND RKGS. 
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EXTERNAL RFX 

REAL Y¥(10) ,DERY (10) ,AC(10) ,v(5) 
COMMON DELT,ERR,BV,V,S,TL,KOU 
INTEGER KOU 

N=NDIM-1 

DO 7 I=1,N 

II=I+1 

DERY (I) -Y (II) 

DERY (NDIM) =RFX (X) 

DO 8 I=1,NDIM 

L=NDIM-I+1 

DERY (NDIM) =DERY (NDIM) -AC (I) *¥ (L) 
RETURN 

END 


SUBROUTINE OUTP (X,Y,DERY,IHLF,NDIM,PRMT,AC, RFX) 
THIS SUBROUTINE COMPUTES THE ISE FROM THE OUTPUT 


VALUES OF RKGS. 


EXTERNAL RFX 


REAL PRMT(6),Y(10),DERY (10) ,AUX(8,10) ,AC(10) ,V(5) 


COMMON DELT,ERR,BV,V,S,TL, KOU 
INTEGER KOU 

IF (IHLF.GE.11)WRITE(6,101) IHLF,X 
FORMAT ('0', 'IHLF=',13,2X,F12.7) 

IF ((PRMT(6)-X).LT.(1.E-5))GO TO 17 
RETURN 

KOU=KOU+1 

ia ROU cEOw os) GOTO 27 

TE (ROUGE 760)GO TO 30 

V (KOU) =RFX (X) 

DO 7 I=1,NDIM,2 

L=NDIM-I+1 

V (KOU) =V (KOU) - (AC (I) *¥ (L) +AC (I) *Y (L) ) 
PRMT (6) =PRMT (6) +DELT 

RETURN 

V (KOU ) =RFX (X) 

DO 8 I=1,NDIM,2 

L=NDIM-I+1 

V (KOU) =V (KOU) - (AC (I) *¥ (L) +AC (I) *¥ (L) ) 
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V1l=V (5) **2 
ERR=ERR+ (BV+3.875* (V(1) **2+V (4) **2)4+2.625* (V(2) **2 
+V(3)**2)+V1) 
IF ( (PRMT (2) -PRMT (6) ) .LT. (1.E-6))GO TO 29 
PRMT (6) =PRMT (6)+DELT 
IF ( (PRMT (2) -PRMT (6)).GE.S)GO TO 31 
KOU=5 
BV=V1 
ERR=ERR* (DELT/3. ) 
RETURN 

81) KOUE0 
BV=V1 
RETURN 

29 ERR=ERR* (DELT/3. ) 
PRMT (5)=0. 
RETURN 

30 Vi=RFX (xX) 

15 DO 9 I=1,NDIM, 2 

L=NDIM-I+1 
9 V1=V1- (AC (I) *Y(L)+AC (I) *Y (L) ) 

Vil=VlLe*2 
ERR=ERR+ (DELT/2. ) * (BV+V1) 
IF ( (PRMT (2) -PRMT(6)).LT.(1.E-6))GO TO 32 
PRMT (6) =PRMT (6) +DELT 
BV=V1L 
RETURN 

32 PRMT(5)=0. 
RETURN 
END 


FUNCTION RFX (X) 
C THIS FUNCTION COMPUTES THE VALUE OF TIME REVERSED SIN 
C (X)/X FUNCTION OF DELAY (2*PI)AND DURACTION (4*PI) 

DOUBLE PRECISION T,TD,TU,P1I,FFX 

TD=2 DO. PL 

TU=4,DO*PI 

Pi=3), 1415927D0 

T=TU-X 

RFX=FFX (T,TD,TU) 

RETURN 

END 


FUNCTION FFX(T,TD,TU) 
C THIS FUNCTION COMPUTS THE VALUE OF SIN(T)/T FUNCTION 
C FOR DELAY (TD) AND DURATION (TU) 

DOUBLE PRECISION T,TD,PI,TX,TU, FFX 

IF(T.GE.TU) GO TO 3 

PI=3.1415927D0 


es 
a eh csyyimesa lt pussys 
eg ot Aeitein el ae 
a 
il 


mat ne 
Os C2) TMS 


ae 
14a) ee eh iar eevee 
' (uvejerty Assia scan 


a as)» “ae haat £) . Wd. 2a) Te (5) 
7 Tws+ ta): 


~ 


BIB-GeeAuySA Sut? dees mv} cone 


(Is) BOLTQARUE | “ 
aa 8 rie wnat ce 


164. 


TX=T-TD 

IE(DABS(TRNELT, V.D-7)) Go TO 1 
FFX=DSIN (TX) /TX 

GO TO 2 

FFX=1.D0 

FFX=FFX/PI 

RETURN 

FFX=0.DO 

RETURN 

END 
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SUBROUTINE ZAGMIN (N,A,E,ALFA,TALF,EXS,KOUNT,A0,A2, 
P,Y,DERY ,AUX) 

THIS SUBROUTINE MINIMISES A FUNCTION OF N VARIABLES 

WITHOUT COMPUTING THE DERIVATIVES.THIS USES THE ZANGWILL 

ALGORITHM. 

SUBROUTINE ZAGMIN(N,A,E,ALFA,TALF,EXS,KOUNT,A0,A2,P,Y, 

DERY, AUX) 

N-#OF VARIABLES-INPUT 

A-VECTOR OF DIMENSION N,CONTAINING THE INITIAL VALUES 

OF THE VARIABLES. 

ON RETURN A CONTAINS THE FINAL VALUES. 

E-FUNCTION VALUE 

ALFA-INITIAL INCREMENT OF SEARCH, INPUT. 

TALF-A SMALL VALUE FOR TESTING THE MINIMUM.IF FUNCTION 

DOES NOT CHANGE FOR ANY DISTANCE LARGER THAN TALF THE 

POINT IS MINIMUM, 

EXS-A SMALL VALUE OF FUNCTION TO TEST MINIMUM.IF 

FUNCTION CHANGE IS LESS THAN EXS IN ONE ITERATION 

MINIMUM IS FOUND 

KOUNT-MAXIMUM # OF ITERATIONS. 

A0O,A2,Y,DERY - EACH OF DIMENSION N IS WORKSPACE 

P-ARRAY OF DIMENSION (N,N) WORKSPACE 

AUX-ARRAY OF (8,N).WORK SPACE 

THIS SUBROUTINE WILL WRITE THE VECTOR A AND FUNCTION 

VALUE AFTER EACH ONE DIMENSIONAL SEARCH. 
RE ATA (0) 2x0 (10) ,A2 (10), ¥ (l0) ,DERY (10) ;AUX (8,10), 
V(5) 
REAL P(10,10) 
COMMON DELT,ERR,BV,V,S,TL, KOU 

INITIALISATION 
N1=N+1 
N2=N-1 
SALF=ALFA 
£C=0 
KO=0 
DOs! Ole =i N 
DO: 204, JI=1,N 
Pinion) ) =O 
Tribe nO. do) P (Ll, Jo)=L. 

101 CONTINUE 
DOMLO2NK= 1 
AO (K) =A (K) 

102 A2(K)=A(K) 
CALL ERROR(N,A, VALUE, Y,DERY , AUX) 
WRITE (6,500)N,ERR, (A(IZAG) , 1ZAG=1,N) 

500 FORMAT ('0',’ORDER',12,1X,'INITIAL ISH',F12.7,1X, 
PAs! ee BFW 2% 7))¥) 
E2=ERR 
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OER=ERR 
C MINIMISATION ALONG CONJUGATE DIRECTIONS. 
1 J=0 
I=0 
2 J=J+1 
IF(J.EQ.N1) GO TO 4 
GO TO 5 
4 J=N 
I=1 
5 DO 103 K=1,N 
103 A(K)=A2 (K)+ALFA*P (K,d) 
WRITE (6,900)J,ALFA 
CALL ERROR(N,A,VALUE,Y,DERY,AUX) 
E3=ERR 
IF (ERR.LT.E2) GO TO 11 
DO 104 K=1,N 
104 A(K)=A2 (K) -ALFA*P (K,J) 
CALL ERROR(N,A, VALUE, Y,DERY,AUX) 
E1=ERR 
IF(ERR.LT.E2) GO TO 8 
IF (ALFA.LT.TALF) GO TO 10 
ALFA=ALFA/2. 
GO TO 5 
8 ALFA=ALFA 
GO TO 11 
10 DO 105 K=1,N 
105 A(K)=A2 (K) 
ERR=E2 
GO TO 15 
11 El=E2 
E2=ERR 
ALFA=ALFA+ALFA 
DO 106 K=1,N 
106 A(K)=A(K)+ALFA*P (K,J) 
WRITE (6,800) ALFA 
CALL ERROR(N,A, VALUE, Y, DERY, AUX) 
TF (PRRILTJE2) "GO HO 11 
E3=ERR 
ALFA=ALFA/2. 
DO 107 OKA LN 
107 A(K)=A(K)-ALFA*P (K,J) 
CALL ERROR(N,A, VALUE, Y,DERY, AUX) 
ir (ERRGUT Pee ico ero TZ 
E3=ERR 
DOs O07 ake, N 
707 A2(K)=A(K) -ALFA*P (K,J) 
DFA= (ALFA/2.) UC (=S2)4+BI-A(4) ). taw-E3yy (bls (2) t2— 
B3:)") 
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ALFA=DFA-ALFA-ALFA 
DO 108 K=1,N 
108 A(K)=A(K)+ALFA*P (K,J) 
CALL ERROR(N,A,VALUE,Y,DERY ,AUX) 
IF(ERR.LE.E2) GO TO 14 
DO 109 K=1,N 
109 A(K)=A2 (K) 
ERR=E2 
GO TO 15 
14 DO 110 K=1,N 
110 A2(K)=A(K) 
E2=ERR 
GO TO 15 
12 E1=E2 
E2=ERR 
DFA= (ALFA/2.)* (((-3.) *El+ (4. ) *E2-E3) / (-El1+ (2. ) *E2- 
E3) ) 
ALFA=DFA-ALFA 
DO 111 K=1,N 
111 A2(K)=A(K)+ALFA*P (K,J) 
CALL ERROR (N,A2,VALUE,Y,DERY, AUX) 
TF(ERR<LT.E2) (GO TO 13 
BORIL2 (KS1,N 
112 A2(K)=A(K) 
ERR=E2 
GO TO 15 
13 E2=ERR 
DOULIsEK= lj N 
113 A(K)=A2 (K) 
15 IF (QeEOUN) BGO TO 16 
19 I=0 
ALFA=SALF 
WRITE (6,300)J3,ERR, (A(IZAG) , 1ZAG=1,N) 
S00 BRORMATUNO', COND DIR. ',l2,1X,!'I1Sb',Fl2,7,1x, A=", 
(8ri2i7)) 
GO TO 2 
16 IF(1.EQ.0)GO TO 17 
TEKABSUOER-ERR) ) LE. (EXS)) GO TO 29 
33 KO=KO+1 
WRITE (6,600) KO,ERR, (A(IZAG) , IZAG=1,N) 
600° FORMAT ('O", "ITERATION toy1x) ISEB’, F12.7,1x, A=", 
Cote 7) ) 
IF (KO.GT.KOUNT) GO TO 30 
IC=IC+1 
ALFA=SALF 
Gourerls 
C SEARCH ALONG N CONJUGATE DIRECTIONS ARE OVER.NEW CON 
C JUGATE DIRECTIONS ARE CHOSEN. 
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17 DO 114 D=1,N 
DO eb4, To=1),N2 
JIJ1=JI+1 
114 P(K,JJ)=P(K,JJ1) 
PL=0. 
DO 115 K=1,N 
P (K,N)=A(K) -AO (K) 
115 PL=PL+P (K,N) *P (K,N) 
PL=SQRT (PL) 
IF(PL.LE.1.E-6) GO TO 33 
DO 116 K=1,N 
116 P(K,N)=P(K,N) /PL 
GO TO 19 
MINIMISATION ALONG THE CO-ORDINATE DIRECTIONS. 
18 I=0 
DO 117 K=1,N 
117 AO(K)=A(K) 
20 I=I+l 
TEGLOSEQeN1) TC=1 
AIC=A(IC) 
21 A(IC)=AIC+ALFA 
WRITE (6,900) IC,ALFA 
900 FORMAT('0',15,2X,'ALFA=',F12.7) 
CALL ERROR(N,A,VALUE, Y,DERY , AUX) 
E3=ERR 
IF (ERR.LT.E2) GO TO 24 
A(IC)=AIC-ALFA 
CALL ERROR(N,A,VALUE,Y,DERY, AUX) 
E1=ERR 
TE(ERR.LT.E2) GO TO 23 
IF (ALFA.LE.TALF) GO TO 22 
ALFA=ALFA/2. 
GOwro: 2 
22° ALIG)=AIC 
WRITE (6,400) 1IC,E2, (A(IZAG) , LZAG=1,N) 
Te(ieHO.N)) GOeTO 31 
IC=IC+1l 
ALFA=SALF 
GO TO 20 
23 ALFA=-ALFA 
24 E1=E2 
E2=ERR 
ALFA=ALFA+tALFA 
A(IC)=A(IC)+ALFA 
WRITE (6,800) ALFA 
800 FORMAT('0', 'ALFA=',F12.7) 
CALL ERROR(N,A, VALUE, Y,DERY, AUX) 
IF (ERR.LT.E2) GO TO 24 
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E3=ERR 
ALFA=ALFA/2. 
A(IC)=A(IC)=ALFA 
CALL ERROR(N,A, VALUE, Y, DERY , AUX) 
IF(ERR.LT.E2) GO TO 26 
E3=ERR 
DFA= (ALFA/2.) * ( ( (-3) *El+(4.) *E2-E3) / (-El+ (2.) *E2-E3) ) 
AIC=A (IC) 
A(IC)=A (IC) -ALFA-ALFA+DFA 
CALL ERROR (N,A,VALUE,Y,DERY, AUX) 
IF(ERR.LT.E2) GO TO 25 
A(IC)=AIC-ALFA 
ERR=E2 
GO TO 28 

25 E2=ERR 
GO TO 28 

26 E1=E2 
E2=ERR 
DFA= (ALFA/2.)*(((-3.) *E1+ (4. ) *E2=E3) / (-El+ (2. ) *E2-E3) ) 
AIC=A (IC) 
A(IC)=A(IC) -ALFA+DFA 
CALL ERROR(N,A,VALUE,Y,DERY,AUX) 
IF (ERR.LT.E2) GO TO 27 
A(IC)=AIC 
ERR=E2 
GO TO 28 

27 E2=ERR 

28 ALFA=SALF 
OER=ERR 
A2 (IC) =A (IC) 
WRITE (6,400) IC,ERR, (A(IZAG) , IZAG=1,N) 

400 FORMAT (10%, 'COOR DIR! ,12,1X,/ISE',F12.7,aXyhA=', 
(emloed )) 
Go TO.1 

C FUNCTION CHANGE IS LESS THAN EXS,ITERATION IS TERMINATED. 

29 WRITE (6,118)EXS 

118 FORMAT('0','ERROR CHANGE IS LESS THAN',F12.7) 
E=ERR 
RETURN 

30 WRITE (6,119) 

119 FORMAT ('0', 'KOUNT IS EXCEEDED') 
E=ERR 
RETURN 

C A MINIMUM IS OBTAINED. 

31 WRITE (6,120) 

120 FORMAT('0','CONVERGENCE IS OBTAINED') 
E=E2 
RETURN 
END 
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SUBROUTINE SERIES (N,A,B,TL,F,X,NT) 
THIS SUBROUTINE COMPUTES THE COEFFICIENTS OF THE 
SERIES EXPANSION OF THE INVERSE OF A RATIONAL LAPLACE 
TRANSFORM 
N-DEGREE OF DENOMINATOR 
A-VECTOR OF DIMENSION N.THIS CONTAINS THE DENOMINATOR 
COEFFICIENTS OF LAPLACE TRANSFORM.A(1)-COEFFICIENT OF 
S**(N-1) ETC. 
B-COEFFICIENT VECTOR OF NUMERATOR,SIMILAR TO A, 
F-A VECTOR TO STORE THE COEFFICIENTS OF SERIES EXPAN 
SION.AN ARBITRARY DIMENSION OF F(300) IS USED.TL- 
MAXIMUM VALUE OF TIME AT WHICH INVERSE IS REQUIRED. 
NT-NUMBER OF TERMS OF SERIES EXPANSION 

DIMENSION A(10),B(10),X(10) ,F (300) 

DOUBLE PRECISION A,B,F,X,FACT,SUM,TL 

TPA ToebE v1. DO) TL=1 D0 

FACT=TL 

NN=N-1 

NNN=N-2 

DO 1 I=1,NNN 

X (I)=0.DO 

X (NN) =TL 

X(N) =-A (1) *TL 

F(1)=B(1) 

DO 6 152,300 

SUM=0.DO 

DO 2 J=1,N 

JI=N-J+1 

SUM=SUM+X (J) *B (JJ) 

F (I) =SUM 

IF((I.GT.N) .AND. (DABS (F(I)).LE.1.D-16))GO TO 7 

FACT=TL/I 

Doug = 1, 

X(J)=X (J) *FACT 

SUM=0. DO 

DO 4 J=1,N 

JJ=N-J+1 

SUM=SUM~-A (JJ) *X (J) 

Dols U-LNN 

JJ=T+1 

X(J)=X (JJ) 

X (N) =SUM 

CONTINUE 

NT=I 

RETURN 

END 
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FUNCTION VAL(NT,F,TL,T) 
C THIS FUNCTION COMPUTES THE INVERSE OF A LAPLACE TRANS 
C FORM USING THE OUTPUT VALUES OF SUBROUTINE SERIES. 
DIMENSION F (300) 
DOUBLE PRECISION F,VAL,X,T,TL 
DEE, DEL DO) TI=1,D0 
X=T/TL 
VAL=F (NT) 
NTT=NT-1 
pO 1 I=1,NTT 
J=NTT-14+1 
VAL=VAL* X+F (J) 
1 CONTINUE 
RETURN 
END 
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COMPUTATION OF THE NRS AND DRS OF THE ORTHONORMAL 
FUNCTIONS FROM THE POLE VALUES,N IS SYSTEM ORDER. IR- 
NO.OF REAL POLES.VECTOR R HAS THE POLE VALUES ,MAGNI 
TUDE OF THE REAL POLES FIRST,FOLLOWED BY COMPLEX POLES, 
MAGNITUDE OF REAL PART FIRST, FOLLOWED BY THAT OF THE 
IMAGINARY PART 
COEFFICIENTS OF NR AND DR OF THE FUNCTIONS ARE PRINTED 
IN ORDER FROM LARGEST TO SMALLEST POWER.ON RETURN R 
VECTOR HAS VALUES AS FOLLOWS:MAGNITUDE OF REAL POLES 
FIRST FOLLOWED BY P1,Q1,P2,Q2 ETC. 
PMPY IS THE SUBROUTINE FOR POLYNOMIAL OPERATIONS GIVEN 
IN SSPLIB OF IBM. 

REAL) Ril0) ,XD( 10) -XNX 20), V3) ,B(10) ,BB(10) /2(10) 

REAL A(10) 

READ(5,1)N 

FORMAT (12) 

READ (5,1)IR 

READ (5,2) (R(T), 1=1,N) 

FORMAT (10F20.7) 

WRITE (6,100) (R(I) , I=1,N) 
GO FORMAT (AO EPOLE VALUES", 10F12,7) 

MOP ys 1 

XN (1)=1. 

TxD=1 

IxNel 

Yi bjy=h. 

y=) 

IF (IR~EQ.0O)GO TO 11 

Vite = 

ty=2 

DOn7ui=171R 

V(hj=R(1) 

Chis DMPYiO5, 02, XD, LXD, ¥,1Y) 

DOga shH las 

Xb (Ry=Z (RK) 

LAD= i Z 

Teviiekoe. tb) GO TO 5 

II=I-1 

Y (1)=-R(II) 

CAL PEPMPY (251 Zy XNPLxNn, Yyly) 

TRN= Z 

DOR4iR=1402 

XN (K) =Z (K) 

DO 6 K=1,1XN 

B(K)=SQRT(R(I)+R(I) ) *XN (K) 

Dou2l KOI=1 ,IXN 

KJ LTI=IXN-KJI+1 

AKI Lb) S]kD (RIIL) 

WRITE (6,101)i, (B(K) ,K=1,1XN) 
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FORMAT (0.7 'NRw OB pd 271K OF 1L 24,7) 
WRITE (6,102) 1, (A(K) ,K=1,1IXN) 
BORMAT (1074 (DROP sy; 124 dx LOW L25:7,) 
CONTINUE 

Y (1)=-R (IR) 

Vi(2)=1, 

Iy=2 

IC=N-IR 

IC=IC/2 

De LO T=1.,0C 

IN=IR+2*I-1 

INN=IN+1 
R(INN)=R(INN) *R(INN) +R(IN) *R(IN) 
R(IN)=R(IN) +R (IN) 

CALL PMPY(Z,1Z,XN,IXN,Y,IY) 
IXN=1Z 

SQRIN=SQRT (2. * (R(IN) ) ) 
SQRINN=SQRT (R(INN) ) 

DO 8 J=1,12Z 

XN (J) =Z (J) 

K=1Z-J+1 

JIJ=JI+1 

B (J) =SQRIN*Z (K) 

BB (JJ) =SQRINN*B (J) 

WRITE (6,203) (XN (KK) ,KK=1,IXN) 
FORMAT ('0','XN',10F12.7) 

IZ=1Z+1 

B(IZ)=0. 

BB(1)=0. 

WRITE (6,103) IN, (B(K) ,K=1,1Z) 
PORMATH. Omy, BOF ,12)1%) 10F 1257) 
WRITE (6,104) INN, (BB(K) ,K=1,12Z) 
EORMAT (201,73 BB OF! 12 ,1X,10F12.7) 
IY=3 

Y (1)=R (INN) 

Y (2)=R (IN) 

¥(3)=1. 

GALL  PMPY (2,,12,XD,,1XD, Y, LY) 
IXD=1Z 

DO 9 k=1,12Z 

XD (K) =Z (K) 

WRITE (6,105) IN, (XD(K) ,K=1,IXD) 
BORMAT(? Ost XD OF 02),.L0P 12.7) 
¥(2)=-¥ (2) 

CONTINUE 

WRITE (6,100) (R(K) ,K=1,N) 

STOP 

END 
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THIS PROGRAMME COMPUTES THE ORTHONORMAL COEFFICIENTS 
AT MINIMUM ISE 
INTEGRATION IS DONE BY USING THE SUBROUTINE DQG32 OF 
SSPLIB. 
FCTX IS A FUNCTION SUBPROGRAMME GIVING THE FUNCTION TO 
BE INTEGRATED 
THIS PROGRAMME WRITES OUT THE ORTHONORMAL COEFFICIENTS 
IN ORDER. 
TL-DURATION OF FUNCTION.THIS SHOULD BE SPECIFIED 
COMMON STATEMENT SHOULD BE USED WITH FUNCTION FCTX 
INPUT DATA IS AS FOLLOWS 
N-ORDER OF ORTHONORMAL FUNCTION 
A-DR OF ORTHONORMAL FUNCTION,A(1) IS THE COEFFICIENT 
OF S** (N-1) 
B-NR COEFFICIENTS SIMILAR TO A 
ORTHONORMAL FUNCTIONS ARE ARRANGED FROM SMALLEST TO 
LARGEST. 
LAST CARD IS BLANK.THIS TERMINATES THE PROGRAMME. 

DIMENSION F(300) ,X(10) ,A(10),B(10) 

DOUBLE PRECISION PI,XL,XU,FCTX,Y,V,T,F,TL,X,A,B 

COMMON F,TL,NT 

EXTERNAL FCTX 

PI=3.1415927D0 

T=PI/2.D0 

TL=4,.DO*PI 

NL=TL/T 

NI=0 

SUM=0.DO 

READ (5,1)N 

IF(N.EQ.0) GO TO 7 

NL=TL/T+.5 

READ (5,2) (A(I) ,I=1,N) 

READS: 2): (B(1) ,I=1,N) 

FORMAT (12) 

FORMAT (10F20.7) 

WRITE (6,3) (A(I) ,I=1,N) 

FORMAT ('0', 'A=',10F12.7) 

WRITE (6,4) (B(I) ,I=1,N) 

FORMAT ('0', 'B=',10F12.7) 

CALL SERIES (N,A,B,TL,F,X,NT) 

XU=0.DO 

V=0.DO 

DO 5 K=1,NL 

XL=XU 

XU=XU+T 

CAnnDOGs2 (Xu, XU, FCTX,Y) 

V=V+¥ 

NI=NI+1 
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WRITE (6,101)NI,V 
FORMAT('0','C',I1,'=',F25.16) 
GO TO. 6 

STOP 

BND 


FUNCTION FCTX (xX) 

DOUBLE PRECISION FCTX,X,PI,F,TL,VAL,FFX 
DIMENSION F (300) 

COMMON F,TL,NT 

PI=3,1415927D0 

TH=2"DO*PI 

FCTX=VAL (NT,F,TL,X) *FFX (X,TD,TL) 

RETURN 

END 


SUBROUTINE MDQG(N,KOUNT,C,EA,GRAD,M) 


IRs) 


THIS SUBROUTINE COMPUTES THE AVERAGE ERROR AND GRADIENTS 
N-ORDER OF SYSTEM. KOUNT-THE ITERATION # TO BE USED BY 
DEMFP.C- COEFFICIENTS OF ORTHONORMAL FUNCTIONS,A N 


VECTOR. EA-AVERAGE ERROR,GRAD-GRADIENTS OF EA WITH 


RESPECT TO THE COEFFICIENTS.M-AN INTEGER,M=0 IF CON 
STRAINT ON ISE IS SATISFIED.M=1 IF THIS IS 


NOT TRUE. 
THIS SUBROUTINE WRITES M,KOUNT,WU,C,EA AND GRAD. 


DOUBLERPRECISION CX(3),X2 (3) FD (10) ,GRAD (10 )\-¢ (10), 


DIMENSION PHI(10,910) 
COMMON WU,WT,RSQ,CH,FD,E EK,XU,T, PHI 
GAMENGOCKNG «CS. XZ) 
NLT= (XU/T)+.5DO 
NLT=NLT*NG* 2 

x= ONO 

EA=0. DO 

IPF=0 

DONe) Ku). 

GRAD (K)=0. DO 

XL=XX 

XX=T+XX 

XA=.5DO* (XX+XL) 
XB=.5DO*T 

DO. 7” K=1) NG 
XT=XZ (K) *XB 

DO 6 KI=1,2 


IPF=IPF+] 

XT=-XT 

X=XA+KT 

R=FFCT (N,C,X,IPF) 

S=SGN (R) 

DO 5 KJ=1,N 

GRAD (KJ) =GRAD (KJ) +CX (K) *FD (KJ) *S 


tC), BRERYXUyT , xX; XL, XA, XB XTX, FECE,5, 2A, SGN, 
*WU ,WT ,RSQ 
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CONTINUE 

EA=EA+CX (K) *DAES (E) 

CONTINUE 

CONTINUE 

IF(IPF.GE.NLT) GO TO 8 

GO TO 4 

EA=XB*EA 

DO 9 K=1,N 

GRAD (K) =XB*GRAD (K) 

CONTINUE 

XxX=0,D0 

DOn12, I=1,N 

XX=XX+ (C(I) -CH(I) ) * (C(I) -CH(I)) 
CONTINUE 

XX=XX-RSQ 

LP Kn UELORDO) mGOnTOugt s 

M=1 

IF (KOUNT.EQ.0) GO TO 14 
EA=EA+WU*XX*XX 

XL=4.DO*WU*XX 

DO 15 -I=1,N 

GRAD (I) =GRAD (I) +XL* (C(I) -CH(I) ) 
CONTINUE 

GO TO 18 

XT=0 DO 

XL=0.DO 

XA=16.DO*XX*XX 

DO 16 I=1,N 
XT=XT+GRAD (I) *GRAD (I) 
XL=XL+XA* (C (I) -CH(I) ) * (C(I) -CH(T) ) 
CONTINUE 

WU=DSQRT (XT/XL) 

XL=4.DO*WU* XX 

EA=EA+WU* XX * XX 

Domb plod, 

GRAD (I) =GRAD (I) +XL* (C (I) -CH(I) ) 
CONTINUE 

WU=WU/WT 

TO TO 18 

M=0 

WRITE (6,10)M,KOUNT ,WU, (C(I) ,T=1,N) 
FORMAT ('0','M=',I1,1X,'KOUNT',14,1X,'WU',F12.4,1X, 
LC=te i SE1247-))) 

WRITE (6,11) EA, (GRAD(I) ,I=1,N) 
FORMAT ('0', 'EA=',F12.7,3X, 'GRAD=', (8F12.7)) 
RETURN 

END 
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SUBROUTINE GQC (NG,CX,XZ) 
C THIS SUBROUTINE GIVES THE CONSTANTS OF 6 POINT GUASS 
C QUADRATURE FORMULA, 
DOUBLE PRECISION CX(3) ,XZ (3) 
NG=3 
XZ (1)=0.93246951420315203D0 
CX (1)=0.17132449237917034D0O 
XZ (2) =0.66120938646626451D0 
CX (2) =0.36076157304813861D0 
XZ (3)=0.23861918608319691D0 
CX (3) =0. 46791393457269105D0 
RETURN 
END 


FUNCTION FFCT (N,C,X,IPF) 
DOUBLE PRECISION C(10),FD(10),CH(10) ,X,E,EK,XU,T, 
*FFCT,FX,WU,WT,RSQ 
DIMENSION PHI(10,910) 
COMMON WU,WT,RSQ,CH,FD,E,EK,XU,T,PHI 
FFCT=0.DO 
DO 1 K=1,N 
FD(K)=PHI(K, IPF) 
FFCT=FFCT+C (K) *FD (K) 
iL CONTINUE 
FFCT=FFCT-FX (X) 
RETURN 
END 


FUNCTION SGN (X) 
DOUBLE PRECISION X,SGN 
PoC ole DO) Orr O = 
Prete Gr .O, 00) GO TO 2 
SGN=0.DO 
RETURN 

i SGN=-1.D0 
RETURN 

2 SGN=1.D0O 
RETURN 
END 


FUNCTION FX (X) 
C THIS FUNCTION GIVES THE VALUE OF DELAYED SIN(X)/X FUNC 
C TION,DELAY=(2*PI), 
C DURATION (4PT) 

DOUBLE PRECISION FX,FFX,P1I,X,TD,TL 

PI=3.1415927D0 

TD=2.DO*PI 

TL=4,DO*PI 
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FX=FFX (X,TD,TL) 
RETURN 
END 


THIS PROGRAMME MINIMISES THE AVERAGE ERROR WITH CON 
STRAINTS ON ISE.THE PENALTY FUNCTION APPROACH IS USED 
TO SATISFY THE CONSTRAINTS.WU IS A WEIGHT FACTOR WHICH 
IS INCREASED BY A FACTOR WT AT THE END OF EACH ITERA 
TION OF DFMFP,IF CONSTRAINT IS NOT SATISFIED. 
E-LEAST ISE AS OBTAINED BY MINIMISATION OF ISE.EK- 
ALLOWABLE ISE FOR MINIMISING THE AVERAGE ERROR.CH- 
COEFFICIENTS OF ORTHONORMAL FUNCTIONS AT LEAST ISE. 
XU-UPPER LIMIT OF INTEGRATION.T-INTERVAL OF INTEGRA 
TION.PHI AN ARRAY OF (10,910) TO STORE VALUES OF ORTHO 
NORMAL FUNCTIONS AT VALUES OF TIME AS REQUIRED BY THE 
GUASS QUADRATURE FORMULA.C-A VECTOR OF DIMENSION N GIV 
ING THE INITIAL VALUES OF COEFFICIENTS FOR MINIMISATION 
OF AVERAGE ERROR. 
THE COMMON STATEMENT IS USED WITH,SUBROUTINES MDQG,FFCT 
AND DFMFP, 
DFMFP IS THE FLETCHER-POWELL MINIMISATION GIVEN IN SSP 
LIB OF IBM. 
A VARIABLE EGS IS ADDED TO THE ARGUMENT LIST OF DFMFP. 
THE MINIMISATION IS TERMINATED WHEN THE SUM OF THE 
MAGNITUDES OF GRADIENTS IS LESS THAN EGS. 
THE STATEMENT IF((M.EQ.1).AND. (KOUNT/N*N.EQ.KOUNT) ) WU= 
WU*WT IS ADDED AT THE BEGINNING OF THE ITERATION OF 
DFMFP 
THIS PROGRAMME SHOULD GIVE THE VALUES OF WT,E,EK,XU,T 
AND INITIAL WU 
THIS PROGRAMME IS SET TO WORK FOR N=T,DELAY=2PI AND 
DURATION=4PI 
INPUT DATA 
N-ORDER OF ORTHONORMAL FUNCTION 
A-DENOMINATOR COEFFICIENTS OF ORTHONORMAL FUNC 
TION.A(1} IS THE COEFFICIENT OF S** (N-1) 
B-NUMERATOR COEFFICIENTS OF THE ORTHONORMAL FUNC 
TIONS.B(1) IS THE COEFFICIENT OF S**(N-1) 
THE FUNCTIONS ARE ARRANGED FROM THE SMALLEST TO LARGEST 
BLANK CARD 
N-ORDER OF APPROXIMATING FILTER. 
CH-ORTHONORMAL COEFFICIENTS AT MINIMUM ISE 
C-INITIAL VALUES OF THE COEFFICIENTS 
DOUBLE PRECISION Cx(3)), x2 (3) ,FD(10),GRAD (L0),C G10), 
*(60),6 (10) ,A(10) ,B(10) ,F (300) ,XD(40)', TL, XV, R, BN, XV, 
aT XX, XL,XA,XB,XT, V, VAL, XU ,WU, WT, RSQ,CH(10) , PI 
DIMENSION PHI(10,910) 
COMMON WU,WT,RSQ,CH,FD,E,EK,XU,T,PHI 
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TL=XU 
PI=3.1415927D0 
EXTERNAL MDQG 
=,2D0 

XU=30.D0 

TL=XU 

CAL ESGOG (NG /CXAXZ) 
NLT= (XU/T) +. 5DO 
NLT=NLT*NG* 2 
IPHI=0 

READ (5,1) N 

FORMAT (12) 
IF(N.EQ.0)GO TO 9 
IPHI=IPHI+1 
READ (5,2) (A(I), 
READ (5,2) (B(I), 
FORMAT (10F20.7) 
IPF=0 

WRITE (6,202) (B(I) ,I=1,N) 
FORMAT ('0', 'B=', 10F12.7) 
WRITE (6,201) (A(I) , I=1,N) 
FORMAT ('0', 'A=',10F12.7) 

CALL SERIES (N,A,B,TL,F,XD,NT) 
X¥X=0.DO 

XL=XX 

XX=XL+T 

XA=. 5DO* (XX+XL) 

XB=.5DO*T 

DO 7 K=1,NG 

XT=XZ (K) *XB 

DO 6 KI=1,2 


I=1,N 


XT=-XT 
XV=XA+t+XT 
IPF=IPF+1 


PHI (IPHI, IPF)=VAL(NT,F,TL,XV) 
CONTINUE 

CONTINUE 
IF(IPF.GE.NLT) GO TO 4 
GO TO 5 

READ (5,1)N 

READ (5,2) (CH(I) ,I=1,N) 
READ (5,2) (C(I) , I=1,N) 
WT=4.DO 

WU=1.D0 

E=,7431E-3 

EK=.035 
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EST=0, 3 

EPS=1.E-4 
LIMIT=100 
EGS=.001 
RSQ=EK-E 
CALL DEMFP (MDQG,N,C,V,G,EST,EPS,EGS, LIMIT, IER,H) 
WRITE (6, 002) (C:( 1)", t= 1) N) 

101 FORMAT('0','C=', (8F12.7)) 
WRITE (6,102) IER,V,G 

VO25 WEORMATE(: Oy abTHR's, LopelX «! BAly bl 2a? tlk, e' GRAD wiGoru 2ue7 ) ) 
STOP 
END 
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103. FORMAT('0','C=',1 
) 
1 


104 FORMAT('0O','R=', 
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THIS PROGRAMME COMPUTES THE NUMERATOR COEFFICIENTS OF 
THE APPROXIMATION 
VECTOR C-COEFFICIENTS OF EXPANSION OF FUNCTIONS 1,2,ETC. 
VECTOR R-MAGNITUDE OF REAL POLES FOLLOWED BY P1,Q1,P2, 
OP Eres 
INPUTS X AND Z ARE THE NR.COEFFICIENTS OF THE RESPEC 
TIVE FUNCTIONS ORDERED FROM SMALLEST TO LARGEST POWER. 
IN DATA NRS ARE READ FROM LARGEST TO SMALLEST FUNCTION 
RESPY. 
INPUT DATA IS ARRANGED AS FOLLOWS. 
N-ORDER OF APPROXIMATION 
IR-NUMBER OF REAL POLES 
C VECTOR AT WHICH NRS OF APPROXIMATION IS DESIRED. 
R VECTOR.MAGNITUDE OF REAL POLES FIRST FOLLOWED BY Pl, 
Oly, P2P oo SIC, 
IX-ORDER OF ORTHONORMAL FUNCTION 
X VECTOR-NRS OF ORTHONORMAL FUNCTION IX ORDERED FROM 
SMALLEST TO LARGEST POWER. 
Z VECTOR-NRS OF ORTHONORMAL FUNCTION (IX-1) ORDERED FROM 
SMALLEST TO LARGEST POWER. THE ORTHONORMAL FUNCTIONS 
ARE READ FROM LARGEST TO SMALLEST RESPY. 
SUBROUTINES PMPY AND PADD ARE FROM SSPUB FOR POLYNOMIAL 
OPERATIONS. 
THE NUMERATOR COEFFICIENTS ARE PRINTED OUT,B(1) COEF 
FICIENT OF S**(N-1) ETC. 
DIMENSION “C(l0),R( 10) © (10), ¥ (10) pm. (10) yen (10) 
READ (5,1)N 
READ(5,1)IR 
FORMAT (I2) 
READ (S29 (C(L), 1 
Ream 2) CR (1), 2 
FORMAT (10F20.7) 
WRITE (6,103) (C(I 


WRITE (6,104) (R(I 


Y(1)=1. 

TY=1 

BAC) = 0", 

IZZ=1 

NX=N-IR 

NX=NX/2 

i Ke NX 
IC=N-2*K+2 

TCC=ic=)) 

READ (5,1) IX 

Pen (5) (x Cl )4 Tak, 1X) 
READ (5,2) G2 (Ll), I=L,2%) 
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DOBSB T= let xX 

X (L)=C (IC) *X (I) +C (ICC) *Z (I) 
CALE QPMPY(Z -1%,%; LX ,¥ , LY) 
DOp4el=1517 

X(I)=Z (I) 

UX= 72 

CALL PADD(Z,1Z,X,IX,22Z,122) 
DOe om ne 

(BAGO VAL) 

IAA Wy: 


- X(1)=R(IC) 


KuGzo=RCLCC) 

3 )=1Ls 

TxX=3 

CALLEEMPY(Z ,12,X,0X%,V;) LY) 
DOG l=. LZ 

Vall iesvad any 

TY=17 

CONTINUE 

TEC TR. ROG) aco. TO 13 
DOeliek=— LR 

IC=IR-K+l 

READ (5,1) 1X 
READ, 2 )iCx (Ll), is Lx) 

DO 8. t= tar 

x(E) ]C(Le)ex (1) 

CALLS PMPV ig ia X EX, phy) 
DOs9> f= bere 

pa eh Ala) 

ix —2 

GAMEPPADO 2 ack, LX, Soph oo) 
POV) 07 

yA es Cl Wea Cal) 

UZ2=L2 

TPC ICIBO. 1) GO TO 13 

TxX=2 

SCLV=R (TC) 

X(2)=1. 

CADGBPMPY (2) 14,4, lh, Xd) 
Done = lene, 

AGO eV aae) 

1A 

CONTINUE 

DOT t=) 2177 

IK=1IZZ-I+1 

SD yHe7 (LK) 

WRITE (6,102) (X(I) ,I=1,122) 
FORMAT ('0', 'B=', 1OF12.7) 
STOP 

END 
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SUBROUTINE FRES (N,A,B,X,DB, PHASE) 
THIS SUBROUTINE COMPUTES THE MAGNITUDE IN DB AND PHASE 
IN RADIANS OF A RATIONAL APPROXIMATION.N-ORDER OF THE 
FILTER.A-COEFFICIENTS OF THE DR. 
A(1) IS THE COEFFICIENT OF THE TERM S**(N-1).B-NR COEF 
FICIENTS SIMILAR TO A. 
X-INPUT FREQUENCY IN RADIANS.DB-MAGNITUDE IN DECIBELS. 
PHASE-ANGLE IN RADIANS 

DIMENSION A(10) ,B(10) 

DOUBLE PRECISION X,Y,Z,XX,YY,2Z2Z,A,B,DB, PHASE 

XX=X 

YY=0. 

ZZ=0. 

Y=B (N) 

Z=A (N) 

DO 7 K=1,N 

IF (K/2*2.EQ.K)GO TO 6 

IF (K.EQ.N)GO TO 10 

KK=N-K 

YY=YY+B (KK) *XX 

ZZ=ZZ+A (KK) *XX 

XX=XX*X 

GO) TO 7 

KK=N-K 

XX=-XX 

IF (K.EQ.N)GO TO 10 

Y=Y+B (KK) *XX 

Z=Z+A (KK) *XX 

XX=XX*X 

CONTINUE 

IF (N/2*2.EQ.N)GO TO 8 

ZL=ZZ+XX 

GO TO 9 

Z=Z+XX 

PHASE=DATAN2 (YY, Y) 

PHASE=PHASE-DATAN2 (ZZ ,Z) 

IF (PHASE.GT.O.DO) PHASE=PHASE-2.DO*3.1415927DO 

Y=Y*Y+YY*YY 

Z=Z*Z+ZZ¥ZZ 

DB=DSQRT (Y/Z) 

DB=20.DO*DLOG1O (DB) 

RETURN 

END 
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APPENDIX II 


LAPLACE TRANSFORMS OF ORTHONORMAL FUNCTIONS 
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